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ABSTRACT 

Smoothed  Particle  Hydrodynamics  (SPH)  is  a  computational  technique  for  the  numerical 
simulation  of  the  equations  of  fluid  dynamics  without  the  use  of  an  underlying  numerical 
mesh.  Although  originally  developed  for  use  in  astrophysical  gas  dynamics,  SPH  has  recently 
been  applied  to  many  other  areas  of  numerical  fluid  dynamics  and  materials  modelling, 
several  of  which  have  particular  relevance  to  defence  problems  of  interest  to  the  DSTO.  In  this 
report  we  review  the  basics  of  the  method  and  then  describe  a  simple  two-dimensional  SPH 
code  for  the  simulation  of  incompressible  fluid  flow.  The  code  is  then  applied  to  simple 
problems  such  as  a  dam  break,  the  sloshing  of  water  and  wave  breaking  over  ships.  These 
examples  illustrate  both  the  capabilities  of  the  technique  and  the  relative  ease  with  which  the 
method  can  treat  problems  which  have  previously  been  considered  difficult  to  solve  using 
traditional  methods  such  as  finite  difference,  finite  volume  or  finite  element  grid  based 
methods.  Further  applications  of  the  method  are  then  reviewed,  concentrating  in  particular  on 
the  utility  of  the  technique  in  solid  mechanics  modelling,  and  then  current  applications  of 
SPH  within  Maritime  Platforms  Division  are  described. 
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Smoothed  Particle  Hydrodynamics: 
Applications  within  DSTO 


Executive  Summary 

Smoothed  Particle  Hydrodynamics  (SPH)  is  a  computational  technique  for  the 
numerical  simulation  of  the  equations  of  fluid  dynamics  without  the  use  of  an 
underlying  numerical  mesh.  Although  originally  developed  for  use  in  astrophysical 
gas  dynamics,  SPH  has  recently  been  applied  to  many  other  areas  of  numerical  fluid 
dynamics  and  materials  modelling,  several  of  which  have  particular  relevance  to 
defence  problems  of  interest  to  the  DSTO.  In  this  report  we  review  the  basics  of  the 
method  and  then  describe  a  simple  two-dimensional  SPH  code  for  the  simulation  of 
incompressible  fluid  flow.  The  code  is  then  applied  to  simple  problems  such  as  a  dam 
break,  the  sloshing  of  water  and  wave  breaking  over  ships.  These  examples  illustrate 
both  the  capabilities  of  the  technique  and  the  relevant  ease  with  which  the  method  can 
treat  problems  which  have  previously  been  considered  difficult  to  solve  using 
traditional  methods  such  as  finite  difference,  finite  volume  or  finite  element  grid  based 
methods. 

We  then  review  the  applications  of  SPH  in  solid  mechanics  modelling.  Examples  of  this 
type  of  work  which  are  of  interest  to  the  defence  community  include  deformation  due 
to  high-velocity  impact,  the  fracture  and  fragmentation  of  cased  explosives,  the 
formation  and  penetration  of  shaped  charge  jets,  and  debris  cloud  dynamics  due  to 
hypervelocity  impact.  SPH  offers  a  new  method  with  unique  capabilities  for  the 
solution  of  such  problems.  We  first  discuss  the  application  of  SPH  itself  to  these 
problems,  then  a  hybrid  approach  which  combines  the  best  features  of  a  Lagrangian 
finite  element  approach  with  the  SPH  method.  This  combined  approach  offers  a 
significantly  improved  capability  for  the  simulation  of  problems  such  as  ballistic 
impact  on  confined  and  unconfined  ceramic  targets  and  transparent  armour. 

A  summary  is  then  given  of  relevant  areas  of  work  in  Maritime  Platforms  Division 
where  SPH  methods  are  currently  being  used.  These  include  the  relative  motion  of  a 
landing  craft  and  the  mother  ship  in  a  well  dock  scenario,  sloshing  within  hulls, 
underwater  explosion  events,  the  deployment  and  retrieval  of  autonomous  vehicles, 
and  ballistic  impact  on  ceramic  targets. 
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1.  Introduction 


Smoothed  Particle  Hydrodynamics  (SPH)  was  developed  independently  by  Gingold  and 
Monaghan  [20]  and  Leon  Lucy  [45]  in  1977  to  simulate  astrophysical  gas  dynamics 
problems.  The  many  techniques  available  at  that  time  for  the  numerical  solution  of  the 
compressible  fluid  dynamic  equations  were  often  unsuitable  for  astrophysical 
applications.  There  were  several  reasons  for  this;  astrophysical  problems  often  involve 
large  changes  in  spatial,  temporal  and  density  scales  over  many  orders  of  magnitude.  Such 
problems  are  also  often  inherently  asymmetric  and  have  no  definite  boundaries.  A  typical 
example  from  astrophysics  which  illustrates  these  problems  is  the  numerical  simulation  of 
the  fission  of  a  rapidly  rotating  star. 

Traditional  techniques  for  the  numerical  solution  of  hydrodynamic  equations  first  involve 
the  creation  of  a  computational  mesh  which  is  used  to  discretise  the  partial  differential 
equations  describing  the  flow.  The  partial  derivatives  can  be  approximated  by  a  number  of 
different  methods  including  finite  difference,  finite  volume,  or  finite  element  schemes.  The 
computational  mesh  can  be  either  fixed  in  space  and  cover  the  entire  fluid  domain  (the 
Eulerian  method),  or  can  be  fixed  to  the  fluid  and  move  with  the  flow  (the  Lagrangian 
method).  The  Eulerian  method  often  requires  the  construction  of  a  very  fine  mesh  over  the 
whole  flow  domain  because  the  location  of  the  interesting  features  of  the  flow  is  not 
known  a  priori.  Hence  the  method  is  often  computationally  expensive.  Fixed  grid  methods 
can  also  suffer  from  excessive  numerical  diffusion  due  to  the  non-linear  convective  terms 
which  arise  in  the  Eulerian  description  of  the  flow. 

Lagrangian  meshes  overcome  these  problems  by  attaching  the  mesh  points  to  the  fluid 
itself  and  allowing  the  points  to  move  with  the  flow.  The  non-linear  terms  then  no  longer 
appear  and  the  mesh  need  only  be  defined  in  the  regions  of  space  occupied  by  the  fluid.  If 
the  motion  of  the  fluid  becomes  geometrically  complex  however  then  the  mesh  undergoes 
severe  distortion  and  the  underlying  numerical  methods  become  unstable  and  the 
computation  stops. 

Neither  of  these  approaches  is  particularly  suited  to  astrophysical  fluid  problems  due  to 
the  lack  of  symmetry,  complex  rotational  behaviour  and  range  of  spatial  scales  inherent  in 
many  of  these  problems.  SPH  was  developed  to  provide  a  computational  scheme  for  the 
solution  of  the  fluid  equations  which  did  not  rely  on  the  use  of  a  computational  mesh.  This 
is  achieved  by  defining  a  set  of  moving  interpolation  points  which  follow  the  fluid  motion. 
In  this  sense  SPH  is  a  Lagrangian  method,  although  the  points  are  never  linked  together  to 
create  a  computational  mesh.  Each  of  the  fluid  dynamic  variables  is  expressed  as  an 
integral  interpolant  using  a  smoothing  function  and  the  integral  is  then  approximated  by  a 
summation  over  the  interpolation  points.  By  using  this  approach,  the  derivatives  of  the 
fluid  variables  can  then  be  evaluated  by  calculating  the  derivatives  of  the  smoothing 
function. 

Although  originally  developed  for  use  in  astrophysical  gas  dynamics,  the  many 
advantages  inherent  in  the  absence  of  a  computational  mesh  have  resulted  in  SPH  being 
adapted  and  applied  to  many  other  areas  of  numerical  fluid  dynamics  and  materials 
modelling.  Typical  applications  include,  free  surface  flow,  high  velocity  impact,  material 
fracture,  detonation  physics,  plasma  dynamics,  marine  fluid-structure  interactions,  and 
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geophysical  fluid  flows.  Many  of  these  areas  of  application  are  relevant  to  research 
currently  undertaken  by  scientists  in  both  the  Maritime  Platforms  Division  (MPD)  and  the 
Weapons  System  Division  (WSD)  of  DSTO. 

The  extent  of  development  and  application  of  SPH  to  the  many  diverse  areas  of  fluid  flow 
and  solids  modelling  is  such  that  no  review  of  SPH  can  hope  to  cover  each  area  of 
application  in  detail.  The  purpose  of  this  report  is  to  provide  a  simple  introduction  to  the 
basics  of  SPH,  to  review  in  some  detail  those  areas  of  application  of  most  relevance  to 
MPD  and  WSD,  and  to  report  simulation  results  from  selected  two-dimensional 
calculations  which  were  performed  to  emphasize  the  attractiveness  and  simplicity  of  the 
method  when  applied  to  problems  which  have  previously  been  considered  difficult  to 
solve  using  traditional  techniques  such  as  finite  difference,  finite  volume  or  finite  element 
grid  based  methods. 


2.  SPH  -  The  Basics 


In  this  section  we  briefly  outline  the  basis  of  the  SPH  method  and  illustrate  the  derivation 
of  the  fluid  dynamic  equations  in  SPH  form.  A  number  of  review  articles  which  describe 
the  foundations  of  the  method  have  already  been  written.  A  relatively  short  paper  written 
by  Monaghan  in  1988  gives  a  clear  derivation  of  the  SPH  equations  and  describes  their 
application  to  a  wide  variety  of  problems  in  compressible  flow  [51].  A  more  detailed  and 
lengthy  review  of  the  method  in  connection  with  astrophysical  gas  dynamics  appeared 
several  years  later  [53].  More  recently,  Monaghan  has  written  a  very  comprehensive 
review  which  details  the  theory  and  application  of  SPH  since  its  inception  in  1977  [55].  As 
well  as  a  detailed  review  of  the  foundations  of  the  subject  Monaghan  also  describes  recent 
applications  of  SPH  in  areas  such  as  incompressible  fluid  flow,  fluid-structure  interactions, 
turbulence,  heat  conduction,  elasticity  and  fracture.  Benz  has  also  written  several  review 
articles  which  outline  the  basis  of  the  method  and  its  applications  in  astrophysics  [2,3]  and 
brittle  solids  modelling  [4].  The  recent  book  by  Liu  and  Liu  [43]  provides  a  very  extensive 
review  of  the  foundations  of  the  subject  as  well  as  providing  many  illustrations  of  the 
application  of  the  method  to  the  simulation  of  detonation,  underwater  explosion  shocks, 
and  the  hydrodynamics  of  material  strength. 

2.1  Integral  Interpolants 

The  basis  of  SPH  is  the  expression  of  each  of  the  fluid  dynamic  variables  as  an  integral 
interpolant.  This  allows  any  function  to  be  expressed  in  terms  of  its  values  at  a  set  of 
disordered  points.  Consider  a  general  function  A(r)  expressed  in  the  form 

4r)  =  |^(r')<j(r-r')dr'  (1) 

where  cj(r-r')  is  the  Dirac  delta  function  and  the  integral  is  taken  over  the  entire  three- 
dimensional  space.  Equation  (1)  is  exact,  but  not  particularly  useful.  The  basic  idea  of  SPH 
is  to  approximate  the  delta  function  by  a  suitable  continuous  function.  In  SPH  this 
function  is  called  the  kernel,  and  the  choice  of  a  suitable  kernel  is  of  central  importance  to 
the  success  of  the  SPH  method. 
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The  integral  interpolant  (A(r))  of  any  function  A(r)  is  defined  as 


(A(y)}  =  f  A(y')jv(y  -Y',/i)dY' 

(2) 

where  the  function  W  has  the  two  properties: 

fA(Y')W(Y-Y',k)dY'  =  l 

(3) 

limW(Y  -Y',k)  =  My-y') 

h^O 

(4) 

Lucy  [45]  referred  to  the  function  W  as  a  broadening  function  because  it  spreads  the 
influence  of  A  on  r  to  the  surrounding  region.  The  length  scale  h  is  effectively  the  half¬ 
width  of  the  kernel  and  determines  the  amount  of  broadening,  or  the  spatial  extent  over 
which  the  function  W  smoothes  the  variable  A.  This  half- width  is  the  resolution  length 
scale  of  SPH  and  is  referred  to  as  the  smoothing  length.  It  is  equivalent  to  the  width  of  a 
grid-cell  in  finite  difference  methods. 

Consider  a  fluid  with  density  p(r)  defined  by  a  set  of  points  rj  initially  distributed  in  a 
regular  manner  throughout  the  body  of  the  fluid.  At  any  time  the  velocity  and  any  other 
fluid  dynamic  variables  are  also  known  at  these  points.  We  can  imagine  the  fluid  to  be 

divided  into  N  small  volume  elements  with  masses  mi,  m2, . ,  mN,  where  the  "centres"  of 

these  small  volumes  are  located  at  the  rf.  Evidently  p(ri)dr  =  mi,  hence  for  numerical 
simulations  the  integral  interpolant  { ^(r))  can  be  approximated  as  follows: 

{A{r))  =  \A^W{r-Y',h)p{r')AY' 

(A(y))  ^J^m.^WiY-Y^h)  (5) 

Pi 

where  the  summation  index  i  denotes  a  particle  label  and  the  summation  is  over  all  N 
particles.  Particle  i  has  mass  mf,  position  rf,  density  pf  and  velocity  vf.  The  value  of  any 
variable  A  at  r/  is  denoted  by  Af. 

The  particular  advantage  of  the  SPH  formulation  becomes  apparent  when  we  consider  the 
integral  interpolant  expression  for  the  gradient  of  a  function 

<  V^(r) )  =  j  VA{y')w{y  -  Y',h)dY'  (6) 

Integrating  by  parts,  equation  (6)  becomes 

<  V^(r) )  =  I  ^(r')VIF(r  -  r',/j)dr'  +  J  A(y')W(y  -  Y')n^.dS  (7) 

The  first  integral  on  the  right  side  of  equation  (7)  is  taken  over  the  volume  of  the  domain, 
while  the  second  integral  is  over  the  boundary  of  the  domain.  In  most  applications  this 
surface  integral  can  be  neglected  because  either  the  function  or  the  kernel  itself  goes  to 
zero  at  the  boundary  of  the  domain.  Equation  (7)  therefore  becomes 
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<v4r))  =  (8) 

which  in  numerical  terms  becomes 

VA(r)  =  ym.^V}F(r-r.,h)  (9) 

/  A 

Equation  (9)  displays  the  fundamental  advantage  of  the  SPH  method;  that  the  derivatives 
of  any  function  A(r)  can  be  found  by  differentiating  the  kernel,  rather  than  by  using  finite 
difference,  finite  volume  or  finite  element  expressions  calculated  from  a  grid. 

The  error  involved  in  approximating  the  function  A(r)  by  the  summation  given  by 
equation  (5)  has  been  studied  in  detail  by  Monaghan  [49].  It  depends  on  the  degree  of 

disorder  of  the  particles,  but  generally  is  of  0(h^)  or  better.  Monaghan  has  discussed  this 
point  in  more  detail  in  his  most  recent  review  [55],  where  he  notes  that  it  is  not  easy  to 
estimate  the  errors  in  the  SPH  equations  from  first  principles  because  the  particles  get 
disordered  during  motion.  The  degree  of  disorder  is  not  complete  however,  as  would  be 
described  by  a  probability  distribution  proportional  to  the  mass  density  for  example, 
because  the  particles  are  still  constrained  by  the  dynamics.  Consequently,  SPH  simulations 
have  been  found  to  be  much  more  accurate  than  the  interpolation  of  quantities  from 
randomly  disordered  particle  arrays  would  suggest.  It  should  also  be  noted  that  although 
the  summation  in  equation  (5)  is  formally  over  all  the  particles,  only  a  small  number 
actually  contribute  because  W  can  be  chosen  so  that  it  falls  off  rapidly  for  |  r  -  rf  |  >h. 

From  a  mathematical  point  of  view  the  points  rf  with  associated  masses  mf  are  simply 
interpolation  points  from  which  the  properties  of  the  fluid  can  be  calculated.  From  a 
physical  point  of  view  however  the  points  can  be  regarded  as  material  particles  with 
masses  mf  which  can  be  treated  like  any  other  particle  system.  It  is  this  physical 
interpretation  of  the  SPH  equations  which  is  most  universally  adopted  and  allows  SPH  to 
provide  a  conceptually  simple  formulation  in  which  the  inclusion  of  more  complicated 
physics  is  a  relatively  straightforward  procedure. 

2.2  Choice  of  kernel  function 

The  original  calculations  of  Gingold  and  Monaghan  [20]  used  a  Gaussian  kernel,  while 
those  of  Lucy  [45]  used  a  bell-shaped  function  constructed  from  a  third  order  polynomial. 
Many  other  functions  are  possible  however,  and  Liu  and  Liu  [43]  give  a  detailed 
discussion  of  the  main  properties  which  any  kernel  function  should  satisfy,  and  also 
provide  a  comprehensive  list  of  some  of  the  most  frequently  used  functions  in  the  SPH 
literature. 

By  far  the  most  popular  kernel  function  is  the  one  devised  by  Monaghan  and  Lattanzio 
[60]  based  on  cubic  spline  functions: 
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=  0  if  q  >  2  (11) 

where  v  is  the  number  of  dimensions  and  a  is  a  normalisation  constant  which  takes  the 
values  2/3,  10 f7n  and  I/tt  in  one,  two  and  three  dimensions  respectively.  This  function 
closely  resembles  the  original  Gaussian  function  used  by  Gingold  and  Monaghan  [20],  but 
has  the  advantage  of  having  compact  support  (meaning  that  the  interactions  are  exactly 
zero  for  r  >  2h)  and  a  continuous  second  derivative.  The  vanishing  of  the  kernel  for  r  >  2h 
is  particularly  important  because  it  means  that  the  summation  in  equation  (5)  needs  to  be 
carried  out  only  for  particles  within  a  radius  of  2h  of  the  point  of  interest. 

Some  of  the  other  functions  described  in  the  literature  include  both  quartic  and  quintic 
spline  expressions  introduced  by  Morris  [61],  which  more  closely  approximate  a  Gaussian 
shape  and  are  more  stable.  Johnson  et  al.  [36]  used  a  quadratic  smoothing  function  because 
it  was  found  to  overcome  a  compressive  instability  problem  which  occurs  in  high  velocity 
impact  problems.  Liu  and  Liu  [43]  recently  introduced  a  new  kernel  based  on  a  fourth 
order  polynomial  expression  which  gives  excellent  results  in  the  standard  shock  tube 
problem  and  in  two-dimensional  heat  conduction  simulations. 

2.3  Equations  of  motion  in  SPH  form 

The  original  SPH  equations  were  derived  to  simulate  astrophysical  gas  dynamic  problems. 
The  momentum  equation,  in  Lagrangian  form,  neglecting  the  effects  of  viscosity  and 
gravity  therefore  takes  the  form 

dv  1 

—  =  —VP  (12) 

dt  p  ^  ^ 

It  is  assumed  that  the  fluid  has  an  equation  of  state  in  which  the  pressure  P  is  a  function  of 
the  density  p  only.  From  equation  (5)  we  can  then  express  the  pressure  gradient  as 

VP  =  Y,m,^VW(r-r„h)  (13) 

b  Pb 

where  we  now  denote  a  particle  label  by  the  summation  index  h.  Hence  equation  (12)  can 
be  written  in  SPH  form  as 


dv 


a 


dt 


j_ 

Pa 


z 


Pb 


(14) 


where  we  now  employ  the  standard  SPH  practice  of  referring  to  the  interpolation  points  as 
particles  and  equation  (14)  is  interpreted  as  the  force  on  particle  a  due  to  all  other  particles 
h  in  the  region  of  particle  fl.  implies  taking  the  gradient  with  respect  to  the  coordinates 

of  particle  a.  and  we  have  now  introduced  the  standard  notation  =  W(rfl  -  tb,  h). 

To  this  we  need  to  add  an  equation  describing  the  motion  of  the  particles  at  each  time  step 


5 


DSTO-TR-1922 


dr, 

dt 


•  = 


(15) 


From  equation  (5)  it  is  also  obvious  that  the  equation  for  the  density  at  any  point  is  given 
by 


Pa 


(16) 


Whilst  the  above  equations  look  reasonable  and  are  the  forms  used  in  the  original  papers 
by  Gingold  and  Monaghan  [20]  and  Lucy  [45],  they  are  usually  modified  slightly  to 
produce  better  results.  Equation  (14)  for  example  does  not  conserve  linear  and  angular 
momentum  exactly  because  the  force  on  particle  a  due  to  particle  h  is  not  equal  to  the  force 
on  particle  h  due  to  particle  a.  The  forces  can  be  made  symmetrical  by  first  writing  VP/  p  in 
the  form 


VP 

P 


P 

+  — V/2 
P 


(17) 


and  then  applying  the  SPH  interpolation  rules  so  that  the  momentum  equation  becomes 


dt 


(4+4>v,w'„ 

Pb  Pa 


(18) 


This  form  of  the  momentum  equation  conserves  linear  and  angular  momentum  exactly.  It 
was  first  derived  by  Gingold  and  Monaghan  using  a  different  approach  [22,23,24]. 

Two  further  modifications  are  particularly  relevant  when  SPH  is  used  to  simulate  the 
motion  of  fluids,  which  is  the  case  in  the  simulations  described  in  this  report.  Instead  of 
moving  the  particles  according  to  equation  (15),  their  motion  is  described  by 


dt 


f , 


=  V 


m. 


Pab  y 


ab 


(19) 


where  =  {pa  +  ph) /  2,  vha  =  -  v^,  and  6*  is  a  constant  with  a  value  normally  of  Vi. 

Equation  (19)  is  known  as  the  XSPH  variant  and  was  introduced  by  Monaghan  [52]  to 
prevent  penetration  when  particle  methods  are  used  to  simulate  streams  of  fluid 
impinging  on  each  other.  The  additional  term  in  equation  (19)  ensures  that  a  particle 
moves  with  a  velocity  that  is  close  to  the  average  velocity  in  its  neighbourhood.  For  the 
simulation  of  nearly  incompressible  fluids  such  as  water  it  has  been  found  to  keep  the 
particles  orderly  in  the  absence  of  viscosity. 

For  many  SPH  simulations  the  density  is  calculated  from  equation  (16).  For  nearly 
incompressible  flow  with  a  free  surface  however,  such  as  water  sloshing  in  a  tank,  use  of 
equation  (16)  produces  incorrect  results.  In  these  types  of  problems  the  density  falls 
discontinuously  to  zero  at  the  surface.  If  equation  (16)  is  used  to  find  the  density  then  the 
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particles  near  the  surface  will  have  their  densities  smoothed  over  a  length  2h,  hence  the 
density  will  not  drop  discontinuously  as  it  should.  The  resulting  smoothing  of  the  density 
discontinuity  produces  incorrect  pressures  and  hence  poor  results  near  the  surface,  which 
in  some  applications  is  the  main  feature  of  interest.  To  overcome  this  problem,  the 
continuity  equation,  i.e. 


^  +  V  •  (/7v)  =  0 
dt  ^  ’ 

can  be  expressed  in  SPH  form  as  follows 

=  \  -V  W 

J ,  "‘a  ^  ab  ^  ab 

dt  V 


(20) 


(21) 


By  using  equation  (21)  instead  of  equation  (16)  the  initial  density  can  be  set  and  then  it  will 
only  change  when  particles  move  relative  to  each  other.  There  is  also  a  computational 
advantage  to  using  equation  (21)  in  place  of  equation  (16)  since  the  rate  of  change  of  all 
physical  variables  can  then  be  computed  in  one  computational  loop. 


2.4  Inclusion  of  viscosity 


Although  Lucy  [45]  was  the  first  to  introduce  a  viscosity  term  in  the  SPH  equations,  the 
most  widely  used  expression  for  viscosity  in  SPH  simulations  was  introduced  by 
Monaghan  and  Gingold  [57].  It  should  be  noted  that  any  numerical  simulation  of  the 
hydrodynamic  equations  which  involve  supersonic  motions  and  the  formation  of  shocks, 
whether  grid  based  or  not,  requires  some  form  of  viscosity  to  stabilize  the  shock.  The 
viscosity  term  introduced  by  Monaghan  and  Gingold  [57]  has  the  advantage  of  containing 
a  term  linear  in  velocity,  which  produces  a  representation  of  the  classical  shear  and  bulk 
viscosities,  as  well  as  a  quadratic  velocity  term  which  handles  high  Mach  number  shocks 
and  is  the  SPH  equivalent  to  the  Von  Neumann-Richtmeyer  artificial  viscosity  used  in 
finite-difference  methods. 


To  understand  the  form  of  the  viscosity  term  as  introduced  by  Monaghan  and  Gingold  we 
first  consider  the  Navier-Stokes  equation  for  viscous  fluids 


dy 

dt 


1  ju  , 

— VP  +  — V"v 
P  P 


(22) 


where  p  is  the  coefficient  of  viscosity.  Note  that  the  addition  of  viscosity  adds  a  second 
derivative  term  to  the  equation.  In  principle  this  is  not  a  problem,  as  second  derivatives 
can  be  estimated  by  differentiating  the  SPH  interpolant  twice.  In  practice  however  the 
expression  can  be  very  sensitive  to  particle  disorder.  To  avoid  this  potential  problem 
Monaghan  and  Gingold  took  a  different  approach.  Equation  (22)  can  be  written  in  the 
following  form 


—  =  V{P-//Vv} 

dt  p  ^  ^  ^ 


(23) 
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which  shows  that  the  viscosity  can  be  regarded  as  an  extra  pressure  term.  Using  this 
analogy,  Monaghan  and  Gingold  took  the  SPH  version  of  the  Navier-Stokes  equation  to 
have  the  form 


^  =  -Z(4  +  4  +  n,.)V,»f^  (24) 

dt  b  Pb  Pa 

where  Ilci,  is  the  additional  "viscous  pressure"  term.  Comparing  equations  (22)  and  (23)  it 
is  obvious  that  IT  should  have  the  form 


n 


(25) 


Now  for  a  gas,  the  coefficient  of  viscosity  ju  is  given  approximately  by 

pXc  (26) 


where  X  is  the  mean-free  path  length  of  the  gas  molecules  and  c  is  the  speed  of  sound  in 
the  gas.  The  logical  extension  of  this  form  of  the  coefficient  of  viscosity  to  SPH  is  simply  to 
replace  X  with  h.  Consider  the  case  of  one-dimensional  motion  first  for  simplicity.  The 
need  to  take  a  second  derivative  of  the  interpolant  can  be  overcome  by  expressing  the 
velocity  derivative  in  equation  (25)  by  a  simple  finite  difference  form 


dy  a  - 

dx  x^-x, 


(27) 


Combining  equations  (25),  (26)  and  (27)  the  expression  for  IT^  (in  one-dimension) 
becomes 


n 


ab 


^  ahc^  ^ 

> 

1 

\  Pa  ) 

(28) 


where  cir  is  a  dimensionless  coefficient  which  can  be  used  to  fine  tune  the  expression  for  a 
particular  simulation.  To  conserve  momentum  however  IT^  needs  to  be  symmetric,  so  Ca 
and  Pa  are  replaced  by  the  symmetrised  forms  and  p^f^ ,  where  )  /  2  and 

Pab  =  {Pa  +  A  2 .  To  avoid  a  problem  when  |  v^b  |  ^  0  and  Xab  approaches  zero  Vab/ Xab 
is  replaced  as  follows 


^ab  _  ^ab  ^ab 


X 


ab 


(29) 


where  rf-  =  0.001/z^.  Generalizing  to  three  dimensions  the  expression  for  IT^  then  becomes 
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n 


ab 


V  P  ab  J 

2  2 

W+V  J 

(30) 


Equation  (30)  successfully  reproduces  the  linear  shear  and  bulk  viscosity  of  fluids  and  was 
also  found  to  stabilize  shocks  of  moderate  strength  [57].  For  astrophysical  problems 
involving  colliding  gas  clouds  however  problems  still  persisted  and  so  an  extra  term  was 
added.  The  additional  term  is  quadratic  in  velocity  and  is  the  SPH  equivalent  of  the  Von 
Neumann-Richtmeyer  artificial  viscosity  used  in  finite-difference  methods.  The  complete 
expression  for  the  general  viscosity  term  IT^  then  takes  the  form 


n 


ab 


Pab 


(31) 


where 


Pab  ~ 


h^ab  .^ab 


ab 


+  7 


When  simulating  high  Mach  number  compressible  gas  flows  where  the  real  viscosity  is 
negligible  then  equation  (31)  only  applies  when  the  fluid  is  being  compressed,  i.e.  when 
particles  are  approaching  one  another,  which  is  the  condition  ®  <  0.  When  ®  > 

0  the  gas  is  expanding  and  then  IT^  =  0.  In  many  applications  to  high  speed  flow  it  has 
been  found  that  the  choice  oi  a=l  and  p=2  produces  excellent  results. 

In  the  opposite  extreme  of  low  Mach  number  flows  where  the  fluid  has  significant  real 
viscosity  then  equation  (31)  can  be  used  with  p=  0.  The  value  of  a  is  then  chosen  so  that 
the  simulated  fluid  has  a  viscosity  and  Reynolds  number  similar  to  that  of  the  real  fluid 
being  modelled.  By  performing  a  more  careful  analysis  of  equation  (31)  Monaghan  and 
Kos  [58]  have  shown  that  it  leads  to  a  shear  viscosity  of  the  form 


ju  =  —ahcp 


(32) 


Monaghan  [54]  has  used  equation  (31)  with  p  =  0  and  a  =  0.01  to  simulate  the  evolution  of 
an  elliptical  water  drop,  the  bursting  of  a  dam,  and  the  formation  of  a  bore  and  in  each 
case  has  found  excellent  agreement  with  experimental  results.  Monaghan  and  Kos  [58] 
also  used  equation  (31)  with  p  =  0  and  a  =  0.001  to  simulate  the  run-up  and  return  of  a 
solitary  wave  travelling  over  shallowing  water  and  then  onto  a  dry  beach  backed  by  a 
vertical  wall.  The  simulations  accurately  reproduced  the  experimental  results.  Equation 
(31)  was  used  to  reproduce  the  effects  of  viscosity  in  all  the  simulations  reported  in  Section 
3  of  this  report. 

Whilst  equation  (31)  has  been  shown  to  give  accurate  results  in  many  scenarios,  including 
typical  fluid  flows  of  interest  in  MPD,  Morris  et  al.  [62]  found  that  it  gave  inaccurate 
velocity  profiles  in  their  simulations  of  Couette  and  Poiseuille  flow  when  the  Reynolds 

number  was  0(10'^).  They  overcame  this  problem  by  using  an  expression  similar  to  the 
first  term  in  equation  (31)  but  based  on  an  SPH  estimation  of  viscous  diffusion  which  is 
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similar  to  an  expression  used  to  model  heat  conduction.  Their  hybrid  expression  still 
combined  a  standard  SPH  first  derivative  with  a  finite  difference  approximation  of  a  first 
derivative.  Takeda  et  al.  [74]  were  the  first  to  use  second-order  derivatives  of  the  kernel  to 
represent  the  second-order  derivatives  in  the  expression  for  the  viscous  term  in  the 
Navier-Stokes  equation.  They  calculated  isothermal  viscous  flow  around  a  cylinder  for 
Reynolds  numbers  in  the  range  6  <  Re  <  55  and  found  excellent  agreement  with  results 
calculated  using  standard  finite  difference  methods.  Watkins  et  al.  [76]  have  shown  how  to 
use  standard  SPH  expressions  for  first  derivatives  to  calculate  accurate  expressions  for  the 
viscous  term  in  the  Navier-Stokes  equation  using  only  first  derivatives  of  the  kernel.  They 
conducted  a  series  of  tests  on  problems  for  which  analytical  solutions  exist  and  found 
excellent  agreement  with  the  theoretical  solutions. 

Balsara  [1]  has  noted  that  the  use  of  equation  (31)  leads  to  substantial  spurious  entropy 
generation  in  regions  of  strong  shear.  He  solved  this  problem  by  multiplying  ]iah  by  a 
suitably  devised  function  of  the  local  compression  and  vorticity  which  approached  unity 
in  regions  of  strong  compression  and  zero  in  regions  of  strong  vorticity.  Calagrossi  and 
Landrini  [11]  slightly  modified  Balsara' s  correction  term  and  used  it  in  a  number  of 
simulations  of  the  classical  dam  break  problem.  They  found  that  this  improved  the  quality 
of  the  solution  both  for  free  surface  and  interfacial  flows,  and  also  improved  total  energy 
conservation. 

For  SPH  simulations  involving  two  or  more  fluids  with  very  different  viscosities  Cleary  [8] 
has  used  a  more  sophisticated  viscous  term  than  that  given  by  equation  (31)  which  allows 
the  viscosity  to  be  variable  and  ensures  that  stress  is  automatically  continuous  across 
material  interfaces.  Cleary  et  al.  [9]  have  shown  that  this  form  of  the  viscosity  allows 
multiple  materials  with  densities  and  viscosities  varying  by  up  to  three  orders  of 
magnitude  to  be  accurately  simulated 

2.5  Treating  incompressible  fluids 

The  development  of  SPH  outlined  so  far  assumes  that  the  fluid  is  compressible.  Since 
water  is  an  almost  incompressible  fluid  some  modifications  are  required  if  SPH  is  to  be 
used  for  hydrodynamic  simulations.  Monaghan  showed  how  this  could  be  achieved  in 

1994  [54].  Since  the  speed  of  sound  in  water  is  of  the  order  of  10^  m/  s  the  Mach  number  in 
typical  simulations  is  extremely  small.  It  is  therefore  customary  in  typical  finite  difference, 
finite  volume  or  finite  element  grid  based  methods  to  approximate  the  fluid  by  an  artificial 
fluid  which  is  exactly  incompressible.  In  SPH  simulations  the  opposite  approach  is  taken; 
the  real  fluid  is  approximated  by  an  artificial  fluid  which  is  more  compressible  than  the 
real  fluid.  The  artificial  fluid  will  still  provide  a  valid  approximation  to  the  motion  of  the 
real  fluid  provided  that  the  speed  of  sound  is  still  much  larger  than  the  speed  of  the  bulk 
flow.  Because  relative  density  fluctuations  are  proportional  to  M^  (where  M  is  the  Mach 
number),  density  fluctuations  can  be  limited  to  -  1%  provided  M  -  0.1.  This  means  that  the 
standard  SPH  formulation  described  above  can  still  be  used.  The  only  modification 
required  is  to  change  the  equation  of  state  of  the  (nearly  incompressible)  fluid  so  that  the 
material  becomes  more  compressible,  while  retaining  sufficient  incompressibility  so  that 
the  speed  of  sound  is  large  enough  to  keep  the  relative  density  fluctuations  small. 
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A  standard  equation  of  state  which  is  used  for  water  in  many  hydrodynamic  simulations 
is  the  one  due  to  Cole  [12],  which  can  be  written  in  the  form 


P  =  B\ 


_P_ 

\Poj 


-1 


(33) 


where  P  and  the  constant  B  are  measured  in  atmospheres  and  po  is  the  density  at 
atmospheric  pressure.  When  n  and  B  take  the  values  7  and  3,000  respectively  this  equation 
agrees  with  the  data  for  water  to  within  a  few  percent  for  pressures  less  than  10^ 
atmospheres.  Monaghan  has  used  this  equation  of  state  to  successfully  simulate  the  flow 
of  elliptical  drops  of  water,  bursting  dams,  the  formation  of  bores,  and  waves  on  a  beach 
[58].  The  only  modification  to  the  equation  of  state  required  was  to  change  the  value  of  the 
constant  B.  With  the  equation  of  state  given  by  equation  (33)  the  speed  of  sound  c  at  the 

reference  density  po  is  given  by 


c 


2 


nB 


Po 


(34) 


Therefore,  if  B  =  lOOpo  v^/n,  the  relative  density  fluctuations  should  be  -  0.01.  This  simply 
requires  that  an  estimate  of  the  maximum  flow  speed  needs  to  be  made  for  each  new 
problem.  A  simple  example  is  a  dam  of  height  H  collapsing.  An  approximate  upper  bound 

on  the  speed  of  the  water  is  then  given  by  v^  =  2gH,  where  g  is  the  strength  of  the 
gravitational  constant. 


2.6  Inclusion  of  boundaries 


Most  boundary  methods  in  SPH  involve  the  use  of  boundary  particles  which  apply  forces 
on  the  fluid  particles.  The  general  form  of  this  force  is  to  treat  the  boundary  force  as  if  it 
were  a  molecular  force.  The  force  is  then  directed  centrally  between  the  particles  and  has 
the  Lennar d-Jones  form 


f(r)  =  D 


(35) 


but  is  set  to  zero  if  r  >  ro  so  that  the  force  is  purely  repulsive.  The  constants  pi  and  p2  must 
satisfy  the  condition  pi  >  p2  and  typically  have  the  values  pi  =  4  and  p2  =  2,  although 

similar  results  are  found  using  pi  =  12  and  p2  =  6.  The  length  scale  ro  is  taken  to  be  the 
initial  spacing  between  the  particles  and  the  coefficient  D  (which  has  dimensions  of 
velocity  squared)  is  again  chosen  by  considering  a  typical  speed  in  the  problem.  For 
example,  for  problems  involving  dams,  bores  or  weirs  with  fluid  of  depth  H  Monaghan 
has  used  D  =  5gH,  but  has  also  shown  that  simulations  using  D  =  lOgH  or  D  =  g¥L  give 
similar  results  [54]. 

Equation  (35)  was  used  by  Monaghan  in  his  initial  free  surface  flow  simulations  [54]  and 
favourable  results  were  obtained.  The  method  is  not  ideal  however  as  it  produces  the 
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equivalent  of  a  corrugated  boundary  containing  ripples  on  the  scale  of  the  particle  spacing. 
Monaghan  and  Koss  [58]  developed  a  better  approach  which  uses  an  interpolation 
procedure  so  that  the  forces  from  neighbouring  boundary  particles  produce  a  force  normal 
to  the  boundary.  In  this  method  boundary  particles  are  assigned  both  a  position  and  a 
local  unit  normal  vector  n  that  points  from  the  boundary  into  the  fluid.  The  force  per  unit 
mass  f  on  a  fluid  particle  from  a  boundary  particle  is  computed  using  the  components  of 
their  separation  along  the  normal  (denoted  by  y)  and  along  the  tangent  (denoted  by  x), 
where  the  distances  x  and  y  are  taken  to  be  positive.  The  force  then  takes  the  form 

f  =  R(y)P(x)n  (36) 

where  R(y)  is  designed  to  fall  to  zero  within  a  few  particle  spacings  of  the  wall.  In  terms  of 
the  variable  q,  where  q  =  y/ (2Ap)  (where  Ap  is  the  initial  particle  spacing),  R(y)  is  defined 

by 

Riy)  =  iiq<l  (37) 

else  R(t/)  =  zero.  The  parameter  A  in  equation  (37)  is  given  by  the  expression 

^  =  |(0.01c" .nj  (38) 

n 


where  the  fluid  particle  is  a  and  the  boundary  particle  is  b.  ^  is  1  if  the  particles  are 
approaching;  otherwise  it  is  zero.  The  second  term  in  equation  (38)  helps  damp-out  the 
motion  perpendicular  to  the  boundary.  The  function  P(x)  is  defined  by 


P(x)  = 


1  +  cos 


V 


if  X  <  Ap 


(39) 


otherwise  P(x)  =  zero.  The  function  P(x)  ensures  that  as  a  fluid  particle  moves  between  two 
boundary  particles  the  contribution  from  the  particles  combines  to  make  the  boundary 
force  constant  if  the  fluid  particle  moves  parallel  to  the  boundary. 

More  recently  Monaghan,  Kos  and  Issa  [59]  have  used  different  expressions  for  the 
functions  R(y)  and  P(x).  If  q  is  now  defined  as  =  y/h,  the  function  R(y)  is  defined  by 


R{y)  =  -P 


V  2  y 

=\pb-qY 


if  0<  q  <  2/3 
if  2/3  <  q  <  1 

if  1  <  q  <  2  (40) 


and  is  zero  for  q  >  2.  The  constant  P  =  O.Olc^f  y  and  is  an  estimate  of  the  maximum  force 
per  unit  mass  necessary  to  stop  a  particle  moving  at  the  estimated  maximum  speed.  The 
factor  1/y  ensures  that  a  particle  moving  faster  than  this  can  also  be  stopped.  This  new 
form  for  R(y)  was  chosen  because  the  boundary  force  opposes  the  pressure  gradient  and 
R(y)  has  the  form  of  the  gradient  of  the  kernel.  The  function  P(x)  is  now  simply  defined  by 
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P(X)  = 


X 


Ap 


if  0<x<Ap 


(41) 


and  is  zero  otherwise.  Monaghan  et  al.  [59]  have  used  the  form  of  the  boundary  force 
defined  by  equations  (40)  and  (41)  in  a  study  of  the  impact  between  a  rigid  body  and  water 
but  have  found  that  there  is  little  difference  between  the  results  obtained  using  these 
expressions  and  the  results  obtained  from  equations  (36)  and  (38).  In  Section  3  we  describe 
the  results  obtained  from  SPH  simulations  using  both  versions  of  the  boundary  force  and 
show  that  the  results  obtained  are  independent  of  the  assumed  form  of  this  force. 


Other  approaches  to  the  problem  of  specifying  realistic  boundary  conditions  are  possible. 
Both  Takeda  et  al.  [74]  and  Morris  et  al.  [62]  performed  simulations  for  low  Reynolds 
numbers  where  the  no-slip  condition  was  the  appropriate  boundary  condition.  To  achieve 
this  they  used  "imaginary  particles"  outside  the  boundary  surface.  Libersky  and  Petschek 
[41]  have  used  a  similar  approach  where  the  additional  particles  are  referred  to  as  "ghost 
particles".  Randles  and  Libersky  [69]  extended  this  approach  to  more  general  boundary 
conditions.  These  schemes  work  by  regarding  the  boundary  as  a  symmetrically  reflecting 
interface,  so  that  the  imaginary  particles  have  the  same  density  and  pressure  as  the 
corresponding  real  particles,  but  opposite  velocity. 


Liu  and  Liu  [43]  use  two  types  of  virtual  particles  to  treat  solid  boundary  conditions. 
Virtual  particles  of  the  first  type  are  located  right  on  the  solid  boundary  and  are  similar  to 
those  used  by  Monaghan  and  described  above.  Virtual  particles  of  the  second  type  are 
similar  to  those  used  by  Libersky  and  Petschek  [41].  Liu  and  Liu  found  that  both  particle 
types  were  needed;  the  virtual  particles  of  the  first  type  ensured  that  the  real  particles 
were  prevented  from  penetrating  the  boundaries,  while  the  virtual  particles  of  the  second 
type  were  needed  to  impose  the  correct  boundary  conditions.  The  various  schemes  differ 
from  one  another  in  the  exact  methods  used  to  calculate  the  velocities  of  the  image 
particles,  the  degree  to  which  these  particles  are  carried  along  by  the  flow  and  are 
included  in  the  summation  process,  and  by  the  way  in  which  the  kernel  sum  is 
normalized. 

In  the  simulations  presented  in  Section  3  we  have  not  used  ghost  particles  nor  any  form  of 
virtual  or  imaginary  particles  outside  the  boundary  surface.  We  have  used  boundary 
particles  as  defined  by  Monaghan  and  Koss  [58]  and  have  included  the  boundary  particles 
in  the  density  calculation  by  solving  the  continuity  equation  in  the  form  of  equation  (21). 
We  note  that  if  fluids  of  very  different  densities  were  present  in  the  simulation  then  a 
different  form  of  the  continuity  equation  would  be  required. 


2.7  Time  stepping  considerations 

The  SPH  formulation  of  the  equations  of  fluid  dynamics  reduces  them  to  a  set  of  ordinary 
differential  equations  for  the  motion  of  each  of  the  particles  within  the  simulation.  Hence 
any  numerical  technique  for  the  solution  of  coupled  ordinary  differential  equations  can  be 
used  for  their  solution.  In  practice  however,  the  right  hand  side  of  the  momentum 
equation  for  each  particle,  equation  (14)  or  equation  (18),  is  generally  quite  expensive  in 
terms  of  computer  time  to  evaluate  and  this  tends  to  exclude  schemes  such  as  high  order 
Runge-Kutta  methods  which  require  several  evaluations  of  the  force  term  at  each  time 
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step.  Most  SPH  codes  use  either  an  improved  Euler  method  (a  mid-point  predictor- 
corrector  method)  [50]  or  a  leapfrog  predictor-corrector  algorithm  for  time  stepping  [51]. 
Each  of  these  methods  requires  only  one  evaluation  of  the  force  term  per  time  step. 

In  the  code  used  to  generate  the  results  shown  in  the  next  section  we  used  the  predictor- 
corrector  leapfrog  algorithm  for  time  stepping.  If  we  write  the  set  of  equations  describing 
the  change  in  velocity,  position  and  density  in  the  following  form 


II 

(42) 

_ 

dt 

(43) 

di 

(44) 

and  denote  the  values  of  the  variables  at  the  beginning  of  a  time  step  At  by  vP, 
and  then  the  predictor  step  is  given  by 

0  T-  0  ,,  0 

1  /  r  ya 

v.,  =  v:+AtF; 

(45) 

F,  +^tv”  +^{htYF^ 

(46) 

p^^=pI  +  md: 

(47) 

New  values  of  Fa  and  Da  are  calculated  using  the  predicted  quantities  and  then  corrected 
values  of  Va  and  p^are  calculated  according  to 


V 


a 


=  V. 


ap 


At(r„-r:) 

(48) 

AtK-A") 

(49) 

The  time  step  is  limited  by  the  familiar  CFL  condition,  which  basically  restricts  the 
physical  rate  of  propagation  of  information  to  be  less  than  that  of  the  numerical 
propagation  rate.  In  SPH  terms  this  becomes  At  <  Iz/c.  If  viscosity  is  present  however  this 
leads  to  an  additional  diffusive  limitation  on  the  time  step  and  the  two  effects  are  usually 
combined  in  the  following  expression 


At 


CV 


min 

a 


h 

c  +  0.6(ac  +  p max^  ) 


(50) 
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A  further  limitation  on  the  time  step  applies  if  external  body  forces  are  present.  This 
implies  that  the  time  step  should  be  less  than  At^^,  where  =  (/z/|f|)^^^  and  f  is  the 

external  force  per  unit  mass  on  each  of  the  SPH  particles.  A  suitable  time  step  for  the 
scheme  therefore  has  the  form 


At  =  ^min(A?„,A^^)  (51) 

The  coefficient  in  equation  (51)  can  be  increased  slightly  whilst  still  maintaining  stability 
and  various  possibilities  are  suggested  in  the  literature  [23,  62]. 

Because  of  the  way  in  which  the  symmetrical  nature  of  the  particle-particle  interactions 
were  maintained  in  the  derivation  of  the  basic  SPH  equations  it  is  possible  to  ensure  exact 
linear  and  angular  momentum  conservation  if  either  a  predictor-corrector  or  leapfrog 
method  is  used  when  solving  the  equations.  Monaghan  has  also  noted  [53]  that,  with  a 
correctly  chosen  time  step,  total  energy  is  conserved  to  within  0.5%  over  400  time  steps. 
Both  integration  methods  are  second  order  accurate  in  time  although  only  one  force 
evaluation  per  time  step  is  required. 

2.8  Variable  smoothing  length 

The  smoothing  length  h  represents  the  effective  width  of  the  kernel  and  its  value 
determines  the  number  of  particles  with  which  a  given  particle  interacts.  The  accuracy  of 
an  SPH  simulation  depends  on  having  a  sufficient  number  of  particles  within  the 
smoothing  length  to  ensure  that  the  replacement  of  an  integral  by  a  summation  is  valid. 
The  speed  of  the  computation  however  decreases  as  the  number  of  such  particles 
increases.  The  optimum  number  of  particles  has  been  discussed  by  Morris  [61]  and 
depends  on  the  number  of  dimensions  in  the  problem.  In  one,  two  and  three-dimensional 
problems  these  are  approximately  5,  21  and  57  respectively.  These  numbers  are  based  on 
the  number  of  neighbours  on  a  cubic  lattice  with  a  smoothing  length  of  1.2  times  the 
particle  spacing  and  a  kernel  which  extends  out  to  Ih. 

If  the  fluid  being  modelled  does  not  undergo  substantial  compression  or  refraction  then  a 
constant  value  for  h  is  sufficient.  However,  if  the  density  changes  substantially  during  the 
course  of  a  computation  then  h  should  be  changed  accordingly  to  maintain  sufficient 
resolution.  In  the  original  formulations  the  resolution  was  constant  in  space  but  allowed  to 
evolve  in  time,  h  =  h{t).  A  method  to  allow  spatially  varying  resolution  was  introduced  by 
Benz  [3]  and  this  resulted  in  enhanced  accuracy  and  speed.  In  such  schemes  each  particle 
has  its  own  h  depending  on  the  local  density  in  the  particle's  neighbourhood  and  the  rate 
of  change  of  the  particle's  density.  For  the  simulation  of  incompressible  fluids  such  as 
water  for  example,  which  are  the  main  interest  of  this  report,  the  density  changes  are 
restricted  to  be  of  the  order  of  1%  and  so  a  constant  value  of  h  was  used  in  all  the 
simulations  described  in  the  next  section. 
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2.9  SPH  coding  details 

There  are  a  number  of  sites  on  the  internet  which  provide  access  to  free  SPH  codes.  These 
include  Joe  Monaghan's  home  page  at  Monash  University^,  which  includes  a  number  of 
sample  FORTRAN  codes  for  the  simulation  of  both  astrophysical  and  fluid  dynamics 
problems,  as  well  as  a  site  at  the  National  University  of  Singapore^,  which  contains  a 
FORTRAN  code  which  was  used  to  perform  many  of  the  simulation  results  discussed  in 
the  book  by  Liu  and  Liu  [43]. 

Any  SPH  code  basically  solves  a  many-body  problem  in  which,  in  principle,  every  particle 
interacts  with  every  other  particle  in  the  problem.  In  practice,  by  using  a  kernel  with 
compact  support,  such  as  that  described  by  equation  (11),  each  particle  interacts  with  only 
a  relatively  small  number  of  neighbouring  particles  confined  within  a  radius  of  2h.  Finding 
the  nearest  neighbours  of  any  given  particle  is  a  major  part  of  any  SPH  code.  The  simplest 
method  of  doing  this  is  to  calculate,  for  a  given  particle  z,  the  distance  rjj  from  z  to  each 
particle  j,  where  j  =  1,2,  ...,  N.  If  the  distance  rjj  is  smaller  than  2h  then  particle  j  is  one  of 
the  nearest  neighbours  of  particle  z.  The  problem  with  this  very  simple  all-pair  search 
approach  is  that  the  time  taken  to  perform  this  search  is  of  order  N^,  and  for  problems 
with  more  than  a  few  thousand  particles  the  computational  time  taken  is  simply  too 
excessive. 

Most  SPH  codes  overcome  this  problem  using  either  a  tree  search  algorithm  or  a  linked  list 
algorithm.  The  tree  search  algorithm  is  popular  in  astrophysical  SPH  codes,  particularly 
those  including  the  effects  of  self-gravity,  and  works  particularly  well  in  problems  in 
which  h  varies  either  in  space  or  time.  This  reduces  the  computational  time  to  0(  N  log  N  ). 
The  linked  list  algorithm  works  well  for  simulations  in  which  h  is  spatially  constant. 
Monaghan  [53]  has  described  the  procedure  for  carrying  out  a  nearest  neighbour  particle 
search  using  linked  lists  in  the  context  of  an  SPH  code.  More  details  can  be  found  in  the 
paper  by  Riffert  et  al.  [70]  and  the  book  by  Hockney  and  Eastwood  [30].  In  essence,  an 
auxiliary  spatial  grid  is  used  to  sort  the  particles  into  cells  and  then  restrict  the  search  to 
neighbouring  cells.  For  example,  if  the  kernel  has  compact  support  of  length  2h  then  the 
mesh  spacing  should  be  set  to  2h.  Then  for  a  given  particle  z  the  nearest  neighbouring 
particles  can  only  be  in  the  same  grid  cell  or  the  immediately  adjoining  cells.  Hence,  in 
three  dimensions,  the  search  is  confined  to  a  maximum  of  27  cells.  The  linked  list 
algorithm  allows  each  particle  to  be  assigned  to  a  cell  and  then  all  the  particles  in  a  given 
cell  are  linked  together  for  easy  access.  The  computational  cost  of  this  algorithm  is 
approximately  O  (N  x  Nneigh)/  where  Nneigh  is  the  average  number  of  contributing 
neighbours  per  particle. 

Another  technique  used  in  SPH  codes  to  reduce  the  computational  time  is  to  use  the 
symmetry  (or  antisymmetry)  inherent  in  the  particle-particle  interactions.  Given  that  the 
force  on  particle  a  due  to  particle  b  is  equal  in  magnitude  (but  opposite  in  direction)  to  the 
force  on  particle  b  due  to  particle  zz,  the  time  to  perform  summations  such  as  those  shown 
in  equation  (18)  can  be  halved.  The  technique  is  to  change  the  range  of  the  indices  in 
nested  DO  loops.  The  symmetry  of  the  interactions  allows  nested  loops  where  both  a  and  b 


^  http ://www.maths .monash. edu. auAjjm/T eaching/welcome.html 
^  http://www.nus.edu.sg/ACES/software/SPH%20code/sph_code_in_the_sph_book.htm 
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go  from  1  to  N  to  be  replaced  by  summations  where  a  goes  from  1  to  N  while  h  goes  from  a 
+  1  to  N.  Additional  computational  time  can  be  saved  by  calculating  the  value  of  the 
kernel  and  its  derivate  over  an  appropriate  spatial  range  at  the  start  of  the  computation 
and  then  storing  the  values  in  a  look-up  table. 


3.  Illustrative  examples  for  fluid  simulations 

In  this  section  we  illustrate  the  basic  SPH  formulation  outlined  above  by  applying  it  to 
some  two-dimensional  problems  of  interest  to  Maritime  Platforms  Division.  The 
simulations  were  performed  using  software  obtained  from  Joe  Monaghan's  website  at 
Monash  University  and  modified  slightly  by  ourselves.  All  simulations  were  performed  on 
a  Pentium  IV  processor  with  640  MB  RAM  running  at  1.80  GHz.  All  programs  are  written 
in  Fortran  and  a  Compaq  Visual  Fortran  6  compiler  was  used  to  run  the  software.  Typical 
run-times  were  no  more  than  30  minutes  using  up  to  20,000  particles. 

3.1  Classical  dam  break 

Monaghan  [54]  used  a  simplified  bursting  dam  problem  as  one  of  several  examples  to 
illustrate  the  ability  of  SPH  to  accurately  model  free  surface  flows.  He  performed  two- 
dimensional  simulations  of  the  collapse  of  a  liquid  column  and  compared  the  results  with 
the  experiments  of  Martin  and  Moyce  [46].  Both  the  height  of  the  dam  and  the  length  of 
the  surge  front  as  a  function  of  time  were  found  to  agree  well  with  experiment,  although 
the  simulated  position  of  the  surge  front  was  found  to  be  slightly  ahead  of  the 
experimental  results  in  all  the  simulations.  Monaghan  attributed  this  to  the  effect  of  drag 
between  the  fluid  and  the  bottom  boundary  particles.  The  boundary  was  simulated  using 
a  line  of  boundary  particles  and  the  force  between  the  boundary  particles  and  the  fluid 
particles  had  the  Lennard-Jones  form  given  by  equation  (35).  Monaghan  noted  that  when 
the  Lennard-Jones  force  was  replaced  by  a  gaussian  force  similar  results  were  obtained. 

Doring  et  al.  [16]  used  both  SPH  and  the  Volume  of  Fluid  (VOF)  method  of  Nichols  and 
Hirt  [63]  to  simulate  the  experiments  of  Martin  and  Moyce  [46]  and  found  similar  results; 
both  the  SPH  and  VOF  simulations  accurately  reproduced  the  height  of  the  dam  as  a 
function  of  time  but  the  surge  front  was  again  slightly  ahead  of  the  experimental  results. 
They  suggested  that  this  difference  was  probably  due  to  the  wall  slip  condition  used  in  the 
calculations.  Colagrossi  and  Landrini  [11]  used  SPH,  a  boundary-element  method,  and  a 
Level  Set  method  to  simulate  the  same  experiments  and  again  found  similar  results;  all  the 
numerical  solutions  agreed  very  well,  but  the  surge  front  was  again  slightly  ahead  of  the 
experimental  results.  They  attributed  the  disagreement  to  experimental  uncertainties  in 
the  early  part  of  the  experiment  and  to  the  neglect  of  important  physical  effects  on  the 
longer  time  scale.  They  suggested  that  an  important  effect  which  was  not  included  in  their 
modelling  was  the  bottom-induced  drag  which  altered  the  propagation  velocity  and 
triggered  the  development  of  turbulence  near  the  water  front.  Liu  and  Liu  [43]  also 
simulated  these  experiments  using  the  SPH  code  described  in  their  book  and  obtained 
results  which  are  more  accurate  than  those  calculated  by  Monaghan  [54]. 
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step  =  200  time  =  0.031 


Figure  1:  Particle  configuration  for  the  collapsing  water  column  at  t  =  0.031  a  =  0.001  and  total 
number  of  particles  =  15,440.  Lengths  are  measured  in  metres 
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Figure  2:  Particle  configuration  for  the  collapsing  water  column  at  t  =  0.506  a  =  0.001  and  total 
number  of  particles  =  15,440.  Lengths  are  measured  in  metres 
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Figure  3:  Particle  configuration  for  the  collapsing  water  column  at  t  =  0.822  a  =  0.001  and  total 
number  of  particles  =  15,440.  Lengths  are  measured  in  metres 
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They  claim  that  this  is  due  to  their  use  of  type  II  boundary  virtual  particles  in  addition  to 
the  type  I  boundary  particles  used  by  Monaghan;  their  use  of  these  additional  virtual 
boundary  particles  in  the  summation  process  increases  the  drag  force  for  the  particles  near 
the  bottom. 

We  have  simulated  the  experiments  of  Martin  and  Moyce  [46]  using  a  two-dimensional 
implementation  of  the  SPH  formulation  outlined  in  Section  2.  Figures  1,  2  and  3  show  the 
particle  configuration  for  the  collapsing  water  column  at  representative  times.  Boundary 
particles  were  used  to  form  the  left-hand  boundary,  the  base,  and  the  right-hand 
boundary.  Three  different  types  of  boundary  force  expressions  were  tried.  Each  of  these 
had  the  form  given  by  equation  (36)  in  Section  2.6  with  the  function  P(x)  given  by  equation 
(39).  The  function  R(y)  had  three  different  forms,  depending  on  the  parameter  potsw.  If 
potsw  =  1  then  R(y)  is  an  integral  repulsive  force  having  the  form  given  by  equation  (37).  If 
potsw  =  2  then  R(y)  is  logarithmically  singular  and  has  the  functional  form  A{l/q  -  1).  If 
potsw  =  3  then  R(y)  is  based  on  the  derivative  of  the  kernel  and  has  the  form  given  by 
equation  (40). 

The  particles  are  initially  set  up  on  a  cartesian  lattice  and  the  smoothing  length  is  defined 
as  /z  =  1.2  Ay,  where  Ay  is  the  initial  particle  spacing.  In  all  the  simulations  described  here  h 
is  constant  in  space  and  time  and  has  the  same  value  for  each  particle.  Due  to  the  manner 
in  which  the  particles  are  initially  placed  in  the  domain  the  system  is  not  in  equilibrium  at 
t  =  0  and  so  a  damping  force  is  applied  to  the  particles  for  the  first  few  hundred  time  steps. 
This  has  the  effect  of  settling  the  system  down  to  an  equilibrium  state  and  is  an  important 
part  of  the  calculation. 


Figures  4  and  5  show  the  simulated  height  and  surge  front  position  as  a  function  of  time. 
All  lengths  are  scaled  by  the  initial  height  of  the  water  column  (HO)  and  time  is  scaled  by 

the  factor  (Ho/g)^/^.  Three  different  runs  are  shown,  corresponding  to  potsw  =  1,  2  and  3. 
In  each  run  a  =  0.001  and  the  total  number  of  particles  was  15,440.  Figure  4  shows 
excellent  agreement  between  the  experimental  results  and  the  simulated  results  for  the 
height  of  the  water  column  as  a  function  of  time.  In  Figure  5  the  agreement  with 
experiment  is  excellent  for  early  times,  but  the  simulated  surge  front  position  leads  the 
experimental  values  at  later  times.  This  behaviour  is  the  same  as  that  found  by  Monaghan 
[54],  Doring  et  al.  [16],  and  Colagrossi  and  Landrini  [11]  and  is  caused  by  the  lack  of  rigour 
in  the  specification  of  the  exact  boundary  condition  at  the  water/ boundary  interface.  The 
effect  is  most  pronounced  at  the  tip  of  the  surge  front  because  the  water  layer  is  very  thin 
at  this  location  and  the  boundary  forces  provide  the  controlling  influence.  Both  Figures  4 
and  5  show  that  the  assumed  form  of  the  boundary  force  has  negligible  effect  on  the 
accuracy  of  the  simulated  results. 
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Figure  4:  Height  of  water  column  as  a  function  of  time,  a  =  0.001  and  the  number  of  particles  = 
15,440 


Scaled  time 


Figure  5: 


Position  of  surge  front  as  a  function  of  time,  a  =  0.001  and  the  total  number  of  particles 
=  15,440 


Figure  6  shows  the  effect  of  changing  both  the  resolution  of  the  calculation  and  the 
assumed  value  of  a  on  the  height  of  the  water  column  as  a  function  of  time.  In  the  figures 
ny  represents  the  number  of  water  particles  initially  placed  next  to  the  left  hand  boundary 
of  the  container.  For  ny  =  80, 120  and  160  the  total  number  of  particles  in  the  calculation  is 
7,120, 15,440  and  26,960  respectively  and  the  corresponding  value  of  is  1.5  cm,  1.0  cm  and 
0.75  cm.  As  can  be  seen,  the  simulated  results  are  virtually  independent  of  the  resolution 
for  the  range  considered  here.  Also,  changing  a  between  a  =  0.001  and  a  =  0.01  has  little 
effect  on  the  calculated  values.  Similar  conclusions  can  be  drawn  from  Figure  7,  which 
shows  the  effect  of  the  above  changes  in  computational  parameters  on  the  simulated  surge 
front  position  as  a  function  of  time. 
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Figure  6:  Height  of  water  column  as  a  function  of  time,  a  =  0.01  and  the  number  of  particles 
varies  between  7,120  and  26,960 
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Figure  7:  Position  of  surge  front  as  a  function  of  time,  a  =  0.01  and  the  total  number  of  particles 
varies  between  7,120  and  26,960 

3.2  Sloshing 

Sloshing  is  defined  as  the  relative  motion  of  a  fluid  with  a  free  surface  confined  inside  a 
tank  caused  by  the  motion  of  the  tank  itself.  It  is  a  highly  non-linear  resonant  phenomenon 
appearing  in  all  marine  structures  containing  liquids  and  is  of  critical  concern  in  the 
design  process  because  its  occurrence  can  lead  to  excessively  high  dynamic  loads  on  the 
tank  walls.  The  simplest  method  to  simulate  the  effects  of  sloshing  is  to  use  an  analogy 
with  a  coupled  mass  and  spring  system,  but  this  does  not  allow  the  pressure  on  the 
container  wall  to  be  calculated.  Solaas  and  Faltinsen  [71]  analysed  the  sloshing  problem 
using  the  potential  flow  hypothesis  by  decomposing  the  free  surface  into  a  sum  of  periodic 
modes.  This  method  works  well  if  the  deformations  of  the  free  surface  are  not  extreme,  but 
is  unable  to  cope  with  events  such  as  wave  breaking  and  roof  impacts,  both  of  which  are 
common  phenomenon  in  ship  tanks. 
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The  availability  of  commercial  finite  element  and  finite  volume  software  packages  such  as 
LS-DYNA,  FIDAP,  Fluent  and  CFX  over  the  past  twenty  years  provides  a  more  accurate 
method  of  simulating  sloshing  phenomena.  Tracking  the  free-surface  of  the  fluid  over  the 
computational  mesh  is  still  a  difficult  problem  however  and  some  of  the  methods  which 
have  been  created  to  overcome  this  problem  include  the  Marker  and  Cell  (MAC)  method 
of  Harlow  and  Welsh  [28],  the  Level  Set  technique  of  Osher  and  Sethian  [67]  and  the 
Volume  of  Fluid  (VOF)  method  due  to  Nichols  and  Hirt  [63].  Whilst  each  of  these  methods 
works  well  in  particular  situations  they  are  computationally  involved  and  expensive.  SPH 
offers  a  much  more  natural  way  of  dealing  with  free-surface  problems  and  it  is  natural  to 
turn  to  this  method  when  attempting  to  simulate  sloshing  phenomena. 

Delorme  et  al.  [14]  used  a  two-dimensional  SPH  method  to  calculate  the  sloshing  pressure 
for  LNG  tankers  and  found  good  agreement  with  results  calculated  using  a  MAC  code. 
Souto  Iglesias  et  al.  [72]  have  used  SPH  to  simulate  sloshing  in  passive  roll-damper  tanks 
for  fishing  vessels.  They  compared  their  results  against  experimental  data  and  found  good 
agreement,  both  for  quantitative  physical  magnitudes  such  as  moment  phase  lags,  as  well 
as  more  qualitative  ones  such  as  free  surface  shapes.  Doring  et  al.  [16]  have  performed 
two-dimensional  SPH  simulations  of  sloshing  in  a  test  tank  and  compared  the  computed 
free  surface  shapes  with  the  experiments  of  Corrignan  [13].  Excellent  agreement  is 
obtained  between  the  simulated  and  experimental  surface  profiles,  even  though 
considerable  non-linearity  is  observed. 

In  this  section  we  have  modified  the  SPH  code  described  in  Section  2  and  used  in  Section 
3.1  to  simulate  the  dam  breaking  problem  by  allowing  the  boundary  particles  which  form 
the  liquid  container  to  oscillate  horizontally.  This  allows  us  to  simulate  the  free  surface 
shapes  obtained  in  the  experiments  of  Corrignan  [13].  In  these  experiments  a  rectangular 
tank  40  cm  wide  and  20  cm  high  containing  water  to  a  height  of  12  cm  is  forced  to  oscillate 
horizontally  with  a  time  dependence  given  by 

X  (t)  =  Ao  [  sin  (27rfit)  -  sin  (27rf2t)  ]  (52) 

where  Aq  =  0.75  cm,  fi  =  1.598  Hz  and  f2  =  1.307  Hz.  The  only  changes  to  the  code  required 
to  do  this  are  to  prescribe  the  motion  of  the  boundary  particles  in  both  the  predictor  and 
corrector  parts  of  the  time  integration  routine.  Initially  we  forced  both  the  position  and 
velocity  of  these  particles  according  to  equation  (52),  but  found  that  this  lead  to 
instabilities  in  the  code  and  resulted  in  some  particles  being  forcefully  ejected  from  the 
fluid  domain.  Stable  behaviour  was  obtained  by  specifying  the  velocity  only,  and  then 
calculating  the  position  from  the  equation 

^«=^:+Atv:  (53) 

Equation  (53)  was  only  applied  during  the  predictor  step,  whilst  the  velocity  was 
prescribed  during  both  the  predictor  and  corrector  step.  This  small  change  in  the 
computational  procedure  lead  to  remarkably  stable  behaviour  and  the  code  was  well 
behaved  over  many  tens  of  thousands  of  time  steps.  We  are  uncertain  as  to  why  this 
particular  procedure  works  as  well  as  it  does,  although  we  note  that  it  is  similar  to  the  type 
of  instability  which  occurs  in  many  CFD  codes  when  certain  boundary  conditions  are  over 
prescribed. 
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It  is  interesting  to  note  that  Gomez-Gesteira  et  al.  [26]  found  similar  behaviour  in  their 
two-dimensional  SPH  code.  They  used  SPH  to  model  wave  overtopping  on  the  decks  of 
offshore  platforms  and  ships  and  used  moving  boundary  particles  to  create  the  initial 
waves.  They  externally  imposed  both  position  and  velocity  on  the  boundary  particles  and 
found  that  this  gave  rise  to  instabilities  and  very  high  instantaneous  accelerations  and 
forces.  They  solved  this  problem  by  imposing  a  smoothing  function  on  both  the  prescribed 
position  and  velocity. 

Figures  8  and  9  show  computed  free  surface  shapes  at  selected  times  for  a  calculation 
using  15,440  particles.  Figure  10  shows  the  comparison  between  the  computed  surface 
shapes  and  the  experimentally  measured  profiles.  The  agreement  is  remarkably  good, 
especially  considering  the  relative  simplicity  of  the  code  employed.  A  similar  calculation 
using  either  a  finite  difference,  finite  volume  or  finite  element  code  would  involve  far 
more  computational  complexity  and  hence  significantly  increased  computational  time.  The 
increased  computational  effort  when  using  a  grid-based  code  results  from  both  the  need  to 
use  specialised  algorithms  to  locate  or  track  the  interface,  as  well  as  the  extra  coding 
involved  to  keep  track  of  the  moving  mesh.  The  amount  of  coding  required  to  implement 
either  of  these  features  in  a  grid-based  code  is  often  greater  than  that  used  to  solve  the 
basic  fluid  equations.  In  an  SPH  simulation,  however,  no  extra  coding  is  required  to  locate 
the  interface  and  in  the  two-dimensional  calculation  shown  here  only  a  few  extra  lines  of 
coding  were  needed  to  effectively  simulate  the  moving  boundaries. 

The  accuracy  of  the  comparison  between  the  computed  and  experimental  surface  shapes 
shown  in  Figure  10  is  typical  of  the  accuracy  of  SPH  simulations.  Other  examples  of  this 
type  of  simulation  are  easily  found.  Next  Limit  Technologies  provide  several  examples  on 
their  website^.  These  simulations  are  conducted  using  their  particle  simulation  code  called 
XFLOW,  which  is  based  on  the  SPH  method,  and  show  two-dimensional  simulations  of 
sloshing  flow  and  comparisons  with  the  experimental  results  by  overlaying  the  simulated 
results  on  the  real  video  footage  of  the  experiment.  The  simulations  run  for  extended 
periods  of  time  and  show  no  degradation  in  the  accuracy  of  the  simulated  surface  shapes 
as  the  length  of  the  simulation  time  increases. 

3.3  Wave  breaking  over  ships 

Water  impact  loading  on  offshore  structures  is  a  subject  area  which  is  now  becoming 
amenable  to  detailed  study  using  sophisticated  computational  fluid  dynamics  codes.  Both 
Neilson  [64]  and  Kleefsen  [38]  have  recently  studied  the  green  water  problem,  where  large 
masses  of  water  can  invade  a  ship's  deck  in  rough  seas  and  cause  considerable  damage  to 
the  ship,  or  to  equipment  or  personnel  on  the  deck.  Both  these  authors  used  finite  volume 
CFD  codes  and  variations  of  the  VOF  method  to  capture  the  free  surface  motion  of  the 
breaking  waves  and  obtained  good  agreement  with  experimental  results,  where  available 
for  comparison. 


^  http://www.nextlimit.com/xflow/index.htm 
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Figure  8:  Simulated  free  surface  shape  for  sloshing  tank  at  t  =  1.94  seconds.  Lengths  are  measured 
in  metres 


step  =  14400  time  =  2.409 


Figure  9:  Simulated  free  surface  shape  for  sloshing  tank  at  t  =  2.409  seconds.  Lengths  are 
measured  in  metres 


Figure  10:  Comparison  between  computed  and  experimental  surface  shapes.  Blue  =  1.65  seconds, 
red  =  2.0  seconds.  The  solid  curves  are  the  SPH  results  while  the  symbols  represent  the 
experimental  results 
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The  use  of  an  Eulerian  grid  to  model  both  the  bulk  of  the  sea  water  and  the  ship  itself 
however  inevitably  results  in  restrictions  on  the  amount  of  the  ship  structure  which  can  be 
included  in  the  simulation.  The  success  of  SPH  demonstrated  so  far  in  accurately 
simulating  free  surface  motion  naturally  leads  to  the  idea  of  combining  SPH  to  model  the 
fluid  motion  with  either  a  finite  volume  or  finite  element  code  to  model  the  ship  structure. 
A  number  of  authors  have  already  addressed  this  problem  and  encouraging  results  are 
being  obtained.  Gomez-Gestiera  et  al.  [26]  have  used  SPH  to  simulate  a  simplified  model 
of  the  green  water  problem.  They  found  that  the  wave  profiles  generated  by  the  method 
were  in  good  quantitative  agreement  with  the  experimental  ones  and  that  the  simulation 
successfully  reproduced  the  main  features  observed  when  a  wave  hits  a  horizontal 
platform.  Gomez-Gesteira  and  Dalrymple  [25]  have  used  a  three-dimensional  SPH  code  to 
model  wave  impact  on  a  tall  structure  and  found  that  the  velocities  and  forces  obtained 
from  their  numerical  code  were  in  excellent  agreement  with  model  laboratory 
measurements.  Fontaine  [19]  has  also  used  SPH  to  model  extreme  waves  and  their 
interaction  with  a  structure  and  Oger  et  al.  [65]  have  used  SPH  to  simulate  wave-body 
interactions  in  extreme  seas.  In  this  section  we  present  a  very  simplified  model  of  a  ship 
sinking  in  rough  seas  to  illustrate  the  capability  of  SPH  to  handle  this  scenario.  The  model 
is  similar  to  one  recently  presented  by  Doring  et  al.  [17].  It  should  be  noted  that  the 
simulation  shown  here  has  been  constructed  to  illustrate  the  capabilities  of  the  SPH 
technique  in  a  simplified  situation;  hence  there  is  no  experimental  data  for  comparison. 

Figure  11  shows  a  sequence  in  which  our  model  ship  (a  rectangular  box  with  a  slit  in  the 
side  to  allow  ingress  of  water)  is  initially  placed  above  the  surface  of  the  water  and  then 
allowed  to  fall  under  gravity.  The  box  then  oscillates  in  a  vertical  plane  until  the  buoyancy 
force  balances  the  gravitational  force  and  the  box  comes  to  a  stable  equilibrium  floating  on 
the  water  surface.  At  time  t  =  2.335  seconds  the  tank  containing  the  water  is  oscillated 
horizontally  to  induce  strong  wave  motion.  Between  t  =  5.0  seconds  and  t  =  5.4  seconds  the 
action  of  the  waves  breaking  over  the  box  results  in  considerable  ingress  of  water  and  the 
box  eventually  sinks  to  the  bottom  of  the  tank  after  7  seconds. 

The  code  used  to  perform  this  simulation  was  the  same  one  which  was  used  to  calculate 
the  results  shown  in  sections  3.1  and  3.2.  To  perform  the  current  simulation  additional 
boundary  particles  were  added  to  create  the  floating  box.  A  subroutine  was  then  written 
which  took  the  total  force  on  the  box  due  to  all  the  fluid  particles  and  then  moved  the 
centre  of  mass  of  the  box  in  accordance  with  this  force.  The  boundary  particles  forming  the 
box  were  then  moved  in  relation  to  the  centre  of  mass  of  the  box  to  simulate  the  motion  of 
a  rigid  body.  The  additional  coding  required  to  perform  this  simulation  was  again 
minimal,  partly  because  we  did  not  allow  the  box  to  rotate  about  its  centre  of  mass.  This 
would  have  involved  some  complicated  geometry  to  calculate  the  new  positions  of  the 
normals  to  the  boundary  particles  and  time  constraints  on  the  work  did  not  allow  this 
level  of  complexity.  Nevertheless,  this  example  clearly  illustrates  the  ability  of  the  SPH 
method  to  cope  with  quite  complicated  fluid-structure  interactions  using  relatively  simple 
coding.  A  comparable  simulation  using  a  grid-based  code  would  involve  an  order  of 
magnitude  increase  in  the  level  of  complexity  of  the  coding  and  the  simulation  time,  for 
reasons  similar  to  those  mentioned  in  the  previous  section.  As  well  as  the  need  to  use 
specialised  algorithms  to  locate  the  interface  and  move  the  mesh  additional  coding  is  also 
required  for  the  fluid-body  interaction. 
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Figure  11:  Sequence  showing  the  sinking  of  an  empty  container  due  to  ingress  of  water  via  wave 
motion.  Lengths  are  measured  in  metres 
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May  and  Monaghan  [47]  have  recently  performed  a  similar  two-dimensional  SPH 
calculation  to  investigate  the  possibility  that  a  large  bubble  of  methane  gas  released  by  the 
ocean  floor  could  capsize  a  floating  body.  In  their  simulation  the  model  ship  was  also 
constructed  from  additional  SPH  boundary  particles  and  was  allowed  to  rotate  as  well  as 
move  in  the  horizontal  and  vertical  directions.  Their  simulations  showed  that  when  the 
radius  of  the  bubble  is  comparable  to  the  length  of  the  ship's  hull  that  it  is  possible  for  the 
bubble  to  cause  the  ship  to  sink.  This  is  due  to  the  mound  of  water  which  is  raised  above 
the  rising  bubble  and  the  subsequent  flow  of  water  from  the  mound  which  creates  a  deep 
trough  either  side  of  the  mound,  which  can  then  carry  the  boat  into  the  trough. 

A  more  sophisticated  approach  to  the  sinking  of  ships  due  to  either  rough  seas  or  bubble 
dynamics  produced  by  underwater  explosions  or  natural  methane  gas  eruptions  would  be 
to  couple  the  SPH  code  to  a  fully  fledged  finite  element  code.  Considerable  progress  along 
these  lines  has  already  been  made  by  Cartwright  et  al.  [5],  [6]  and  [7].  These  authors  have 
used  an  SPH  model  embedded  in  the  commercial  finite  element  code  PAM-SHOCK  to 
simulate  the  non-steady  motion  of  various  vessels  in  varying  sea  conditions.  The 
simulations  to  date  have  used  rigid  finite  element  models  for  the  vessels  so  that  accurate 
force  magnitudes  could  not  be  obtained.  Nevertheless,  the  simulated  vessel  motion 
displays  qualitatively  correct  behaviour  for  all  of  the  vessels  studied.  A  simulation  of  the 
Incat  91  metre  Wave  Piercing  Catamaran  showed  that  dynamic  lift  from  the  under  body 
shape  of  the  hull  was  developed  as  the  correct  cruise  speed  was  approached  [5],  while  in  a 
further  application  the  method  was  shown  to  account  for  the  correct  dynamic  loads  on  a 
34  metre  catamaran  sailing  yacht  [6].  The  most  advanced  application  of  this  method  to 
date  has  been  to  the  simulation  of  the  motion  of  a  landing  craft  within  the  flooded  well 
dock  of  a  parent  ship  [7].  Very  encouraging  preliminary  results  have  been  obtained  and 
this  example  will  be  discussed  in  more  detail  in  Section  5. 

4.  Applications  to  solid  mechanics  modelling 

Problems  involving  large  scale  material  deformation  have  always  been  of  interest  to  the 
defence  community.  Examples  include  deformation  due  to  high-velocity  impact,  the 
fracture  and  fragmentation  of  cased  explosives,  the  formation  and  penetration  of  shaped 
charge  jets,  and  debris  cloud  dynamics  due  to  hypervelocity  impact.  The  simulation  of 
these  problems  has  traditionally  been  conducted  using  large  computer  programs  known 
as  hydrocodes,  which  solve  the  continuum  equations  for  the  conservation  of  mass, 
momentum  and  energy  in  combination  with  sophisticated  equations  of  state  which 
describe  material  behaviour. 

During  the  1980s  and  1990s  WSD  used  a  number  of  these  codes  to  assist  in  the  design  of 
various  warheads  and  the  analysis  of  their  behaviour.  The  Lagrangian  finite  element  code 
DYNA2D  was  used  to  model  the  formation  of  Explosively  Formed  Projectiles  (EFPs)  [32], 
Shaped  Charge  jets  [33],  the  investigation  of  limpet  mine  damage  effects  against  ship  hulls 
[75],  and  the  design  of  an  EFP  to  be  used  as  a  stand-off  sea  mine  neutralisation  device  [40]. 
The  use  of  Lagrangian  codes  to  simulate  material  response  is  limited  to  cases  where  the 
amount  of  deformation  is  small.  When  the  deformation  is  large,  Eulerian  codes  must  be 
used  and  WSD  used  the  finite  volume  Eulerian  codes  HELP  [27]  and  HULL  [18]  in  these 
situations.  These  codes  can  handle  arbitrarily  large  material  deformations  but  require 
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sophisticated  coding  techniques  to  model  material  interfaces.  Difficulties  can  also  occur 
when  simulating  problems  with  fracture,  void  formation,  and  fragmentation. 

During  the  period  1990-1996  defence  scientists  in  both  the  USA  and  UK  became  aware  of 
the  unique  capabilities  of  the  SPH  method  and  began  to  adapt  it  to  the  solution  of  solid 
mechanics  problems  relevant  to  the  defence  community.  Although  the  name  Smoothed 
Particle  Hydrodynamics  might  be  thought  to  imply  that  only  hydrodynamic  problems 
could  be  addressed  by  the  method,  material  strength  can  in  fact  be  added  in  a  fairly 
straight  forward  manner.  Libersky  and  Petschek  [41]  were  the  first  to  show  how  this  could 
be  done  by  incorporating  an  elastic-perfectly-plastic  material  strength  model  into  the  SPH 
framework  and  applying  it  to  the  simulation  of  the  impact  of  an  iron  rod  with  a  rigid 
surface.  Early  results  were  quite  encouraging  and  the  method  was  quickly  extended  into 
the  three-dimensional  shock  and  material  response  code  MAGI  [42],  which  was  based 
entirely  on  the  SPH  method.  The  code  was  used  to  model  additional  cylinder  impact  tests 
and  the  debris  cloud  resulting  from  the  hypervelocity  impact  of  a  copper  disc  on  an 
aluminium  plate.  A  companion  paper  using  two-dimensional  calculations  [68]  (therefore 
achieving  a  significantly  finer  resolution)  of  the  same  hypervelocity  impact  problem 
showed  superb  agreement  between  the  SPH  particle  plot  and  the  X-ray  photographs  of  the 
debris  cloud. 

An  excellent  summary  of  the  early  work  of  Libersky,  Petschek  and  co-workers  can  be 
found  in  the  paper  by  Randles  and  Libersky  [69].  This  paper  also  highlights  the  significant 
advantages  of  the  SPH  method  for  simulating  the  dynamic  response  of  materials  involving 
fracture  and  fragmentation.  Several  simulations  of  the  detonation  of  cased  explosives  are 
given,  including  a  simulation  of  the  detonation  of  a  MK82  general  purpose  bomb  out  to 
280  ps.  The  fragment  distribution  data  was  found  to  compare  remarkably  well  with  the 
experimental  data. 

It  should  be  stressed  that  a  very  important  property  of  the  SPH  method  is  that  it  allows  a 
totally  seamless  transition  between  a  fully  continuum  description  of  matter  and  the 
transformation  of  the  material  into  any  number  of  fragments  produced  by  impact.  At  least 
one  noted  authority  in  this  area  believes  that  it  is  definitely  the  best  method  for  dealing 
with  the  fragmentation  of  brittle  solids  by  impact  [56].  The  only  underlying  problem  with 
the  application  of  SPH  in  this  area  is  that  the  normal  SPH  elastic  equations  do  not  conserve 
angular  momentum.  This  has  been  noted  by  Monaghan  [56]  and  discussed  by  Hoover  et 
al.  [31],  who  have  shown  that  by  using  strong  XSPH  smoothing  the  loss  in  angular 
momentum  of  a  rotating  elastic  wheel  can  be  considerably  reduced. 

A  more  common  approach  to  the  use  of  SPH  methods  in  the  simulation  of  problems 
involving  large  scale  material  deformation  and  fracture  is  to  couple  the  SPH  method  with 
a  standard  Lagrangian  finite  element  code.  The  Lagrange  method  (where  the  mesh  distorts 
with  the  material)  has  the  advantage  of  being  fast  and  providing  excellent  definition  of  the 
material  interfaces.  When  the  deformation  becomes  too  large  however  the  grid  becomes 
entangled  and  the  calculation  stops.  Lagrangian  codes  try  to  overcome  this  problem  by 
using  an  erosion  algorithm,  where  elements  which  have  reached  a  specified  value  of  strain 
(typically  150%)  are  simply  removed  from  the  grid.  This  is  obviously  a  non-physical 
process  however  and  results  in  energy  being  artificially  removed  from  the  calculation. 

SPH  provides  the  perfect  solution  to  this  problem;  as  the  elements  become  highly  strained 
they  are  simply  converted  to  SPH  nodes.  This  immediately  removes  the  problem  of  grid 
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entanglement  as  the  SPH  particles  are  never  required  to  have  any  connectivity.  Johnson  et 
al.  [34]  were  the  first  to  show  how  this  could  be  achieved  using  the  Lagrangian  finite 
element  code  EPIC.  A  two-dimensional  axisymmetric  version  of  the  code  was  used  and 
simulation  results  were  shown  for  both  cylinder  impact  and  hypervelocity  impacts  with  a 
variety  of  different  material  models  including  copper,  iron,  steel,  aluminium  and  concrete. 
Comparison  with  experimental  results  was  excellent.  More  details  of  the  method,  and 
more  simulation  results  for  high  velocity  impact  problems,  are  provided  in  Johnson  [35] 
and  Johnson  et  al.  [36]. 

Shortly  after  the  work  of  Johnson  appeared  Swegle  and  Attaway  [73]  coupled  the  SPH 
method  into  the  transient  dynamics  finite  element  code  PRONTO  and  used  it  to  simulate  a 
number  of  underwater  explosion  problems  involving  fluid-structure  and  shock-structure 
interactions.  Their  results  showed  that  this  combined  method  was  well-suited  to  model  the 
transmission  of  loads  from  underwater  explosions  to  nearby  structures.  Bubble  formation 
and  collapse  was  also  modelled  effectively,  although  the  authors  noted  that  the  method 
still  had  difficulty  (like  all  previous  methods)  in  calculating  the  late  time  effects  due  to  the 
acceleration  of  gravity  and  bubble  buoyancy. 

Following  on  from  the  work  of  Johnson  et  al.  [34]  and  Swegle  and  Attaway  [73],  Hayhurst 
et  al.  [29]  implemented  an  SPH  module  in  the  two-dimensional  axisymmetric  version  of 
the  AUTODYN  hydrocode  for  the  simulation  of  ballistic  impact  problems.  AUTODYN  is  a 
fully  integrated  suite  of  codes  incorporating  Lagrange,  Shell,  ALE  (Arbitrary  Lagrange 
Eulerian)  and  Eulerian  solution  techniques  which  can  be  coupled  in  several  ways  in  space 
and  time.  Previous  numerical  simulations  of  ballistic  impact  problems  using  AUTODYN 
used  either  the  Lagrange  or  Euler  processors,  both  of  which  have  significant  failings  in 
certain  areas.  As  noted,  the  Lagrangian  method  suffers  from  grid  entanglement  problems 
when  the  material  deformation  becomes  large.  This  problem  can  be  overcome  using  the 
Eulerian  approach,  where  the  grid  remains  fixed  in  space,  but  the  codes  are  considerably 
more  expensive  to  run  due  to  the  complicated  algorithms  required  to  resolve  material 
interfaces,  and  the  method  is  also  ill-suited  to  the  implementation  of  sophisticated  fracture 
mechanics  models,  which  require  the  complete  history  of  the  material  to  be  followed 
(which  is  the  case  with  the  Lagrangian  approach). 

Hayhurst  et  al.  [29]  tested  their  SPH  version  of  AUTODYN-2D  by  simulating  the  impact  of 
an  iron  cylinder  on  a  rigid  wall,  a  steel  projectile  impacting  a  ceramic  target  backed  by 
aluminium,  and  the  penetration  of  a  tungsten  long  rod  into  a  thick  steel  target.  In  each  case 
the  SPH  simulation  performed  as  well  or  better  than  simulations  using  either  the  Euler  or 
Lagrange  processors.  They  concluded  by  noting  that  the  combination  of  an  SPH  algorithm 
with  a  Lagrangian  processor  maintained  all  the  advantages  of  the  Lagrangian  method, 
such  as  the  efficient  tracking  of  material  interfaces  and  the  ability  to  incorporate 
sophisticated  material  models,  while  removing  the  problem  of  grid  tangling  in  a 
physically  realistic  manner. 

Clegg  et  al.  [10]  further  pursued  this  approach  by  using  the  SPH  capability  in  AUTODYN- 
2D  to  simulate  kinetic  energy  penetrator  impacts  on  multi-layered  soil  and  concrete 
targets.  The  advantages  of  the  technique  were  clearly  highlighted  in  this  series  of 
simulations  because  the  extent  of  cracking  in  the  concrete  was  far  more  realistically 
simulated  using  the  combined  Lagrangian/ SPH  method.  Neither  the  Eulerian  nor 
Lagrangian  calculations  were  as  effective  in  predicting  spall  from  both  the  front  and  back 


29 


DSTO-TR-1922 


sides  of  the  target.  The  use  of  the  Lagrangian/SPH  technique  allowed  sophisticated 
constitutive  models  to  describe  the  concrete  behaviour,  such  as  hydrostatic  compaction, 
yielding,  damage  and  cracking,  which  are  extremely  difficult  to  implement  in  Eulerian 
codes,  while  overcoming  the  problem  of  mesh  tangling.  The  method  proved  to  be  so 
successful  that  the  SPH  method  was  implemented  in  AUTODYN-3D  and  was  used  to 
simulate  the  formation  of  an  Explosively  Formed  Penetrator  (EFP),  a  Shaped  Charge  Jet, 
and  hard  penetrator  impact  onto  ceramic  armour  [66].  This  version  of  the  code  has  also 
been  used  to  model  the  impact  of  stainless  steel  and  tantalum  projectiles  onto  transparent 
targets  at  speeds  which  are  sufficiently  high  enough  to  result  in  the  shattering  of  the 
projectile. 

In  the  discussion  regarding  wave  breaking  over  ships  in  Section  3.3  it  was  noted  that  the 
finite  element  code  PAM-SHOCK  contained  an  SPH  module  and  that  the  code  had  been 
used  to  simulate  the  non-steady  motion  of  vessels  in  varying  sea  conditions.  In  that 
application  of  the  code  the  SPH  module  was  used  to  model  the  bulk  motion  of  the  water. 
Kamoulakos  et  al.  [37]  have  used  the  SPH  module  in  PAM-SHOCK  for  a  completely 
different  application.  They  conducted  a  number  of  space  debris  impact  simulations  on 
Whipple  shields  using  both  the  SPH  option  and  the  finite  element  version  of  the  code  and 
found  good  agreement  with  experiment  using  both  methods. 

The  Livermore  Software  Technology  Corporation  software  LS-DYNA  is  another  example 
of  a  commercial  finite  element  code  which  has  recently  been  coupled  with  an  SPH  module. 
Details  of  the  implementation  can  by  found  in  the  paper  by  Lacome  [39],  which  also 
contains  examples  of  the  application  of  the  combined  finite  element/ SPH  solver  to  the 
simulation  of  bird  impacts  on  turbine  blades  and  various  types  of  ballistic  impact. 

De  Vuyst  et  al.  [15]  have  recently  implemented  an  SPH  algorithm  in  the  public  domain 
version  of  the  DYNA  code.  Their  method  of  implementation  uses  a  contact  force  vector  to 
treat  the  finite  element  nodes  as  SPH  particles.  This  is  different  to  the  earlier 
implementations  of  SPH  in  finite  element  codes,  which  typically  used  a  master-slave 
algorithm  to  couple  the  two  techniques.  The  performance  of  their  algorithm  is  illustrated 
by  the  simulation  of  three  diverse  impact  problems:  a  plate  impact,  water  impact  and  rod 
penetration.  Excellent  agreement  was  found  between  the  simulation  results  and  the 
experimental  or  numerical  results  for  each  of  these  problems. 

5.  Current  SPH  work  in  MPD 

Several  research  groups  within  MPD  are  currently  using,  or  planning  to  use,  the  SPH 
technique  for  the  simulation  of  a  number  of  different  problems.  These  include  the  relative 
motion  of  a  landing  craft  and  the  mother  ship  in  a  well  dock  scenario,  sloshing  within 
hulls,  underwater  explosion  events,  the  deployment  and  retrieval  of  autonomous  vehicles, 
and  ballistic  impact  on  ceramic  targets.  This  section  briefly  describes  each  of  these  areas  of 
application. 

5.1  Surface  Platform  Systems  Branch 

The  examples  described  in  Section  3  of  this  report  show  that  SPH  provides  a  convenient, 
simple  and  robust  method  for  the  study  of  free  surface  fluid  motion.  It  is  this  aspect  of 
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SPH  which  makes  it  a  particularly  useful  tool  with  several  areas  of  application  within  the 
Surface  Platform  Systems  (SPS)  Branch. 

One  of  the  current  major  projects  for  the  ADO  is  the  acquisition  of  two  large  amphibious 
Landing  Helicopter  Dock  (LHD)  ships  containing  well  docks.  In  order  for  DSTO  to  be  able 
to  provide  advice  to  minimise  the  risk  associated  with  potential  operational  constraints 
during  the  selection  of  the  final  design  there  is  a  need  to  be  able  to  simulate  the  relative 
motion  between  the  LHD  and  the  landing  craft  within  the  well  dock.  Using  conventional 
finite  element,  finite  difference  or  finite  volume  codes  this  would  be  a  problem  of 
considerable  difficulty  and  limited  application.  The  combination  of  a  finite  element  code  to 
model  the  LHD  and  the  landing  craft  with  an  SPH  module  to  model  the  fluid  motion 
however  makes  this  problem  much  more  amenable  to  simulation. 

Considerable  progress  in  this  area  has  already  been  made  by  Cartwright  and  McGuckin  of 
Pacific  ESI  in  collaboration  with  Stuart  Cannon  and  Terry  Turner  from  the  SPS  branch. 
They  used  the  Pacific  ESI  finite  element  code  PAM-SHOCK,  which  contains  an  embedded 
SPH  module,  to  simulate  the  motion  of  two  generic  landing  craft  models  within  the  well 
dock  of  a  parent  ship  [6].  Figure  12,  an  illustration  of  the  capability  of  the  code,  shows  a 
landing  craft  about  to  enter  the  flooded  well  dock  of  a  parent  ship. 

Simulations  were  run  for  two  scenarios.  In  the  first  the  landing  craft  was  tethered  in  a 
fixed  position,  while  in  the  second  the  landing  craft  was  moving  forward  within  the  well 
dock.  The  results  showed  that  this  method  provided  a  viable  simulation  tool  for  the 
prediction  of  the  relative  motion  between  the  LHD  and  the  landing  craft,  although  it  was 
noted  that  several  research  challenges  needed  to  be  solved  and  experimental  validation 
was  required  before  the  method  could  be  considered  ready  for  application  to  a  specific 
problem. 

One  of  the  problems  with  the  SPH  method,  noted  by  Cartwright  et  al.  in  a  previous  paper 
[6],  is  that  considerable  loss  of  wave  amplitude  can  occur  in  an  SPH  simulation  if  a  wave  is 
propagated  over  distances  of  more  than  several  wavelengths.  This  problem  is  currently 
being  addressed  by  Cannon  and  Turner  in  collaboration  with  Joe  Monaghan  (Monash 
University)  and  Paul  Cleary  (CSIRO).  Experimental  wave  data  obtained  by  Turner  in  the 
towing  tank  at  the  Australian  Maritime  College  in  Launceston  is  being  used  to  benchmark 
simulations  conducted  by  Monaghan  and  Cleary  using  in-house  SPH  codes.  Once  the 
cause  of  the  energy  loss  is  identified  and  appropriate  solution  strategies  are  implemented 
it  is  anticipated  that  these  refined  SPH  methods  will  be  implemented  in  future  SPH-finite 
element  code  combinations  and  lead  to  more  accurate  predictive  capabilities.  Current 
work  in  this  area  recently  reported  by  Monaghan  [56]  has  shown  that  by  using  a  different 
time  stepping  algorithm  there  was  little  decrease  in  amplitude  over  five  wavelengths,  but 
that  wave  heights  had  reduced  to  about  60%  of  the  initial  amplitude  after  propagating 
over  ten  wavelengths.  This  problem  does  not  occur  with  the  simulation  of  solitary  waves 
generated  by  a  numerical  wavemaker. 
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Figure  12:  Simulation  of  the  docking  of  a  landing  craft  within  an  LHD  ship  using  the  FAM- 
SHOCK  finite-element/SPH  software 


Another  problem  of  current  interest  within  the  SPS  branch  concerns  sloshing  within  ship 
hulls.  Interest  in  this  area  has  been  stimulated  by  the  recent  grounding  of  HMS 
Nottingham  on  Wolfe  Rock  near  Norfolk  Island.  The  ingress  of  water  and  the  internal 
sloshing  caused  by  the  motion  of  the  grounded  ship  due  to  the  impact  of  the  waves  lead  to 
significant  weakening  of  the  internal  bulkheads.  Initial  simulation  work  in  this  area  is 
currently  being  done  using  the  LS-DYNA  finite  element  code,  which  also  contains  a 
coupled  SPH  module  for  the  simulation  of  fluid  motion.  As  a  precursor  to  the  simulation 
of  a  complete  hull  DYNA  is  currently  being  used  to  simulate  the  sloshing  experiments  of 
Neilson  [64].  Results  to  date  show  that  the  code  can  accurately  calculate  the  shape  of  the 
fluid  surface,  but  some  problems  are  being  experienced  in  calculating  accurate  pressure 
loadings  on  the  tank  walls. 

This  is  not  an  uncommon  situation  with  current  SPH  calculations.  The  method  works 
remarkably  well  in  areas  remote  from  boundaries,  but  special  effort  is  required  to 
accurately  represent  forces  on  boundaries.  As  mentioned  in  Section  3.3,  Gomez-Gesteira 
and  Dalrymple  [25]  have  used  a  three-dimensional  SPH  code  to  model  wave  impact  on  a 
tall  structure  and  found  that  pressure  loadings  obtained  from  their  numerical  code  were  in 
excellent  agreement  with  model  laboratory  measurements.  They  achieved  this  result  by 
treating  the  boundary  particles  as  fixed  fluid  particles.  The  advantage  of  this  approach  is 
that  it  is  unnecessary  to  make  an  a  priori  assumption  about  the  nature  of  the  force  exerted 
by  the  boundaries,  nor  to  utilize  special  computations  there;  a  repulsive  force  normal  to 
the  boundary  results  without  any  additional  considerations. 

Once  the  boundary  force  problems  in  DYNA  have  been  overcome  it  is  anticipated  that  the 
code  will  then  be  used  to  simulate  loads  on  ship  bulkheads.  This  information  will  then  be 
used  as  input  to  short  term  damage  assessment  codes  currently  in  use  in  MPD.  This  will 
give  an  assessment  of  the  strength  of  the  vessel  in  such  situations  and  aid  in  the  initial 
application  of  damage  control  procedures. 
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As  mentioned  in  Section  4.0,  SPH  has  already  been  applied  to  the  simulation  of 
underwater  shock  and  bubble  dynamics  produced  by  the  detonation  of  underwater 
explosives.  An  understanding  of  these  events  and  an  ability  to  predict  their  behaviour  is 
vital  for  the  prediction  of  warship  survivability.  DSTO  has  therefore  been  engaged  in  an 
experimental  and  computational  study  of  underwater  explosions  for  some  years  now,  with 
the  ultimate  aim  of  developing  a  realistic  simulation  capability  for  their  effects  on  naval 
vessels.  To  date  this  work  has  been  carried  out  with  both  Lagrangian  and  Eulerian  finite 
element  codes.  These  mesh-based  techniques  have  limited  success  in  following  the 
changing  bubble  geometry  over  the  course  of  several  bubble  expansions  and  contractions 
however  and  so  some  preliminary  work  is  being  done  on  this  problem  using  the  SPH 
method.  This  is  a  collaborative  project  between  Dr.  Irene  Penesis  at  the  Australian 
Maritime  College  in  Launceston  and  the  group  led  by  Dr.  John  Brett  in  MPD.  Initial 
studies  will  be  performed  using  the  SPH  code  described  in  the  book  by  Lui  and  Lui  [43], 
suitably  modified  for  the  modelling  of  underwater  explosives.  The  project  will  simulate 
the  motion  of  explosive  bubbles  and  compare  the  simulated  results  with  experimental  data 
obtained  by  the  MPD  group  from  experiments  conducted  in  Epping  quarry. 

5.2  Undersea  Platform  Systems  Branch 

An  ideal  tool  for  maritime  automation  research  has  recently  been  developed  in  the  UPS 
branch  by  John  Wharington.  The  software  is  known  as  Odessa  and  provides  an  integrated 
simulation  environment  for  vehicular  and  articulated  systems.  It  simulates  the  motion  of 
multiple  rigid  bodies  interlinked  by  various  types  of  joints  and  has  applications  in  many 
maritime  scenarios  including  simulation  of  AUV,  ROV  and  UAV  dynamics,  payload 
drops,  decoys  and  towed  arrays.  Odessa  is  based  on  ODE  (Open  Dynamics  Engine),  an 
open-source  package  used  primarily  for  game  development.  The  software  uses  the  LCP 
(Linear  Constraint  Problem)  formulation  but  makes  various  approximations  in  order  to 
produce  faster  simulations  with  greater  robustness. 

Odessa  incorporates  a  special  cable  modelling  package  developed  at  RMIT  University  over 
the  past  decade  which  is  specific  to  DSTO  requirements.  This  specialized  cable  code  is 
faster  and  more  accurate  than  using  primitive  elements  and  is  capable  of  explicit 
incorporation  of  bending,  torsional  stiffness,  and  hydrodynamic  loads,  thus  making  it 
ideal  for  the  simulation  of  umbilical  cables  on  ROVs  and  AUVs,  multipart  tow  systems, 
and  towed  arrays. 

To  extend  the  capabilities  of  Odessa  to  allow  the  code  to  model  complex  multibody 
interactions  involving  fluids,  a  new  fluids  simulation  module  has  been  written  and 
integrated  into  the  software  by  Royce  Smart.  Because  the  applications  envisaged  include 
highly  non-linear  free  surface  fluid  motion  the  software  uses  the  SPH  method.  The  code  is 
based  on  the  SPH  program  contained  in  the  book  by  Liu  and  Liu  [43]  but  has  been 
completely  re-written  in  C++  to  more  easily  interface  with  the  Odessa  code  and  the 
graphics  software.  The  use  of  the  SPH  technique  is  also  appropriate  because  the  method  is 
remarkably  robust  and  can  be  adjusted  to  provide  a  trade  off  between  speed  and  accuracy, 
which  is  also  a  feature  of  the  underlying  LCP  method  in  Odessa. 

The  incorporation  of  the  SPH  module  in  Odessa  will  allow  the  code  to  model  well-dock 
scenarios  and  the  deployment  and  retrieval  of  autonomous  vehicles.  This  work  is 
complementary  to  the  well-dock  simulations  conducted  in  the  SPS  branch  conducted  by 
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Cannon  and  Turner  using  PAM-SHOCK  because  the  sophisticated  cable  model  in  Odessa 
will  allow  more  realistic  simulations  of  the  various  methods  of  restraining  the  landing 
craft  within  the  well-dock,  while  the  finite  element  capabilities  in  PAM-SHOCK  will 
allows  more  realistic  simulations  of  individual  types  of  landing  craft. 

The  current  capability  of  the  SPH  module  within  Odessa  is  illustrated  in  Figure  13,  which 
shows  a  sequence  in  the  damn  break  problem  at  intervals  of  0.5  seconds.  The  code  is  fully 
three-dimensional  and  the  illustrations  show  an  SPH  fluid  interacting  with  Odessa  rigid 
body  objects,  which  are  the  transparent  walls  of  the  container.  The  SPH  particles  are  colour 
coded  by  velocity,  with  blue  being  the  lowest  velocity  and  red  the  highest  velocity. 

5.3  Advanced  Materials  and  Sensors  Branch 

MPD  has  core  responsibility  for  research  into  enhancing  military  vehicle  survivability  and 
has  a  long  history  of  achievement  in  ballistic  protection  research  in  all  areas  of  armour 
technology.  Predictive  modelling  plays  a  significant  role  in  evaluating  the  performance  of 
various  types  of  armour  against  ballistic  impact  and  previous  work  in  this  area  has  been 
conducted  using  the  Lagrangian  finite  element  code  DYNA. 

The  prevalence  of  new  composite  armours  containing  ceramic  layers  however  has  lead  to 
a  need  to  upgrade  current  simulation  capabilities  to  deal  with  the  complicated  fracture 
mechanics  of  these  types  of  materials.  Stephen  Cimpoeru  and  his  team  in  the  AM&S 
branch  are  planning  to  use  the  SPH  module  in  the  AUTODYN  software  package  to 
simulate  the  ballistic  protection  offered  by  confined  and  unconfined  ceramic  targets  and 
transparent  armour.  As  discussed  in  Section  4.0,  the  SPH  module  within  AUTODYN  has 
already  displayed  considerable  success  in  simulations  of  this  type.  The  combination  of  the 
finite  element  Lagrangian  solver  with  the  SPH  module  and  the  sophisticated  materials 
models  within  AUTODYN,  which  includes  constitutive  models  for  metals,  composites, 
ceramics,  glass,  concrete,  soil,  and  explosives,  currently  provides  the  most  appropriate 
tools  for  modelling  ballistic  impact  events  of  this  type. 

6.  Discussion  and  Conclusion 


The  use  of  the  SPH  technique  for  the  simulation  of  hydrodynamic  flows,  gas  dynamics 
and  solids  modelling  has  increased  significantly  during  the  past  decade.  The  areas  in 
which  SPH  has  found  application  are  now  so  diverse  that  no  review  can  authoritatively 
cover  each  of  these  areas  in  significant  detail.  The  simulations  and  applications  discussed 
in  this  report  have  concentrated  on  those  areas  thought  to  be  of  most  use  to  current 
defence  science  interests. 

A  large  body  of  work  has  now  shown  that  SPH  is  a  particularly  useful  tool  for  the 
prediction  of  bulk  fluid  motion  with  free  surfaces.  SPH  is  clearly  the  best  numerical 
technique  to  apply  to  these  types  of  problems  if  the  simulated  results  for  the  fluid  motion 
are  not  required  to  display  extremely  high  levels  of  accuracy,  which  is  typically  the  case  in 
many  defence  applications.  Examples  of  this  type  of  application  include  the  sloshing  of 
water  inside  hulls,  the  well  dock  problem,  and  wave  breaking  over  ships. 
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Figure  13:  A  sequence  from  a  damn  break  problem  at  0.5  second  intervals  calculated  using  the 
Odessa  code  and  the  three-dimensional  SPH  module 
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However,  one  area  of  hydrodynamics  in  which  SPH  does  not  offer  any  specific  advantage, 
and  in  fact  is  definitely  not  the  most  appropriate  solution  technique,  occurs  in  simulations 
where  viscous  forces  are  the  dominant  factors  in  determining  the  nature  of  the  flow. 
Examples  of  this  type  of  application  include  the  calculation  of  lift  and  drag  forces  on 
submerged  bodies  in  high  Reynolds  number  flows.  In  these  applications  the  flow  is  highly 
turbulent  and  the  general  features  of  the  flow  are  determined  by  the  detachment  of  the 
boundary  layer.  The  dynamics  of  this  are  highly  dependent  on  the  viscosity  of  the  fluid. 
SPH  can  model  drag  forces  successfully,  as  the  papers  by  Takeda  et  al.  [74]  and  Morris  et 
al.  [62]  have  shown,  but  only  for  very  low  Reynolds  number  flows  which  are  essentially 
laminar  and  where  the  boundary  layer  can  be  resolved  by  the  SPH  particles.  In 
applications  where  the  Reynolds  number  is  of  the  order  of  a  million  or  higher  the 
thickness  of  the  boundary  layer  is  several  orders  of  magnitude  smaller  than  the  typical 
length  scale.  This  causes  significant  problems  for  any  numerical  solution  technique 
because  of  the  need  to  resolve  length  scales  over  several  orders  of  magnitude.  An  SPH 
code  with  a  variable  smoothing  length  could  certainly  be  applied  to  such  problems,  but 
would  currently  offer  no  specific  advantages  over  traditional  grid  based  methods. 

The  other  area  where  SPH  has  had  a  significant  impact  in  an  area  of  science  relevant  to 
defence  is  in  the  simulation  of  the  brittle  fracture  of  solids.  This  is  because  the 
sophisticated  fracture  mechanics  models  used  to  accurately  simulate  the  behaviour  of 
these  materials  requires  the  entire  stress  history  of  a  given  piece  of  material.  This  is  easily 
done  using  Lagrangian  codes  because  in  this  approach  the  frame  of  reference  is  attached  to 
the  material.  Lagrangian  methods  however  suffer  from  the  problem  of  grid  entanglement 
and  become  impractical  to  use  in  many  applications  of  defence  interest.  The  ability  to 
replace  entangled  Lagrangian  nodes  with  disconnected  SPH  particles,  as  discussed  in 
section  4.0,  means  that  the  Lagrangian  approach  can  now  still  be  applied,  even  in  cases  of 
extreme  material  distortion.  As  noted  in  section  5.3,  several  commercial  codes  now  contain 
this  technique  for  fracture  mechanics  modelling  and  this  offers  a  significantly  improved 
capability  for  the  simulation  of  ballistic  impact  on  confined  and  unconfined  ceramic 
targets  and  transparent  armour. 
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